| blob | d3cafd48d1869d12581e53de21d7a19c1cebeba3 |
1 !WRF:MODEL_LAYER:PHYSICS\r
2 !\r
3 MODULE module_bl_acm\r
4 \r
5 ! USE module_model_constants\r
6 \r
7 REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number\r
8 REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER\r
9 \r
10 CONTAINS\r
11 \r
12 !-----------------------------------------------------------------------\r
13 !-----------------------------------------------------------------------\r
14 SUBROUTINE ACMPBL(XTIME, DTPBL, U3D, V3D, &\r
15 PP3D, DZ8W, TH3D, T3D, &\r
16 QV3D, QC3D, QI3D, RR3D, &\r
17 #if (WRF_CHEM == 1)\r
18 CHEM3D, VD3D, NCHEM, & ! For WRF-Chem\r
19 KDVEL, NDVEL, NUM_VERT_MIX, & ! For WRF-Chem\r
20 #endif\r
21 UST, HFX, QFX, TSK, &\r
22 PSFC, EP1, G, &\r
23 ROVCP, RD, CPD, &\r
24 PBLH, KPBL2D, EXCH_H, EXCH_M, REGIME, &\r
25 GZ1OZ0, WSPD, PSIM, MUT, RMOL, &\r
26 RUBLTEN, RVBLTEN, RTHBLTEN, &\r
27 RQVBLTEN, RQCBLTEN, RQIBLTEN, &\r
28 ids,ide, jds,jde, kds,kde, &\r
29 ims,ime, jms,jme, kms,kme, &\r
30 its,ite, jts,jte, kts,kte)\r
31 !-----------------------------------------------------------------------\r
32 !-----------------------------------------------------------------------\r
33 \r
34 ! THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO \r
35 ! THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2 (ACM2), WHICH IS A COMBINED \r
36 ! LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)\r
37 !\r
38 ! REFERENCES: \r
39 ! Pleim (2007) A combined local and non-local closure model for the atmospheric\r
40 ! boundary layer. Part1: Model description and testing. \r
41 ! JAMC, 46, 1383-1395\r
42 ! Pleim (2007) A combined local and non-local closure model for the atmospheric\r
43 ! boundary layer. Part2: Application and evaluation in a mesoscale\r
44 ! meteorology model. JAMC, 46, 1396-1409\r
45 !\r
46 ! REVISION HISTORY:\r
47 ! AX 3/2005 - developed WRF version based on the MM5 PX LSM\r
48 ! RG and JP 7/2006 - Finished WRF adaptation\r
49 ! JP 12/2011 12/2011 - ACM2 modified so it's not dependent on first layer thickness. \r
50 ! JP 3/2013 - WRFChem version. Mixing of chemical species are added\r
51 ! HF 5/2016 - MPAS version \r
52 ! JP 8/2017 - Z-coord version from MPAS version for hybrid coords \r
53 ! JP 10/2020 - Fixed a bug where KSRC is now an array dimensioned ( its:ite )\r
54 ! - Modified the PBL height algorithm for stable conditions so\r
55 ! that the Richardson number is computed using windspeed in layer\r
56 ! k rather than wind speed difference between layer k and ksrc\r
57 ! - Added EXCH_M to argument list \r
58 !\r
59 !**********************************************************************\r
60 ! ARGUMENT LIST:\r
61 !\r
62 !... Inputs:\r
63 !-- XTIME Time since simulation start (min)\r
64 !-- DTPBL PBL time step\r
65 !-- U3D 3D u-velocity interpolated to theta points (m/s)\r
66 !-- V3D 3D v-velocity interpolated to theta points (m/s)\r
67 !-- PP3D Pressure at half levels (Pa)\r
68 !-- DZ8W dz between full levels (m)\r
69 !-- TH3D Potential Temperature (K)\r
70 !-- T3D Temperature (K)\r
71 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)\r
72 !-- QC3D 3D cloud mixing ratio (Kg/Kg)\r
73 !-- QI3D 3D ice mixing ratio (Kg/Kg)\r
74 !-- RR3D 3D dry air density (kg/m^3)\r
75 !-- CHEM3D Chemical species mixing ratios (ppm) Optional for WRFChem \r
76 !-- VD3D Dry deposition velocity (m/s) Optional for WRFChem\r
77 !-- NCHEM Number of chemical species Optional for WRFChem\r
78 !-- UST Friction Velocity (m/s)\r
79 !-- HFX Upward heat flux at the surface (w/m^2)\r
80 !-- QFX Upward moisture flux at the surface (Kg/m^2/s)\r
81 !-- TSK Surface temperature (K)\r
82 !-- PSFC Pressure at the surface (Pa)\r
83 !-- EP1 Constant for virtual temperature (r_v/r_d-1) (dimensionless)\r
84 !-- G Gravity (m/s^2)\r
85 !-- ROVCP r/cp\r
86 !-- RD gas constant for dry air (j/kg/k)\r
87 !-- CPD heat capacity at constant pressure for dry air (j/kg/k)\r
88 !-- GZ1OZ0 log(z/z0) where z0 is roughness length\r
89 !-- WSPD wind speed at lowest model level (m/s)\r
90 !-- PSIM similarity stability function for momentum\r
91 !-- MUT Total Mu : Psfc - Ptop\r
92 \r
93 !-- ids start index for i in domain\r
94 !-- ide end index for i in domain\r
95 !-- jds start index for j in domain\r
96 !-- jde end index for j in domain\r
97 !-- kds start index for k in domain\r
98 !-- kde end index for k in domain\r
99 !-- ims start index for i in memory\r
100 !-- ime end index for i in memory\r
101 !-- jms start index for j in memory\r
102 !-- jme end index for j in memory\r
103 !-- kms start index for k in memory\r
104 !-- kme end index for k in memory\r
105 !-- jts start index for j in tile\r
106 !-- jte end index for j in tile\r
107 !-- kts start index for k in tile\r
108 !-- kte end index for k in tile\r
109 !\r
110 !... Outputs: \r
111 !-- PBLH PBL height (m)\r
112 !-- KPBL2D K index for PBL layer\r
113 !-- REGIME Flag indicating PBL regime (stable, unstable, etc.)\r
114 !-- RUBLTEN U tendency due to PBL parameterization (m/s^2)\r
115 !-- RVBLTEN V tendency due to PBL parameterization (m/s^2)\r
116 !-- RTHBLTEN Theta tendency due to PBL parameterization (K/s)\r
117 !-- RQVBLTEN Qv tendency due to PBL parameterization (kg/kg/s)\r
118 !-- RQCBLTEN Qc tendency due to PBL parameterization (kg/kg/s)\r
119 !-- RQIBLTEN Qi tendency due to PBL parameterization (kg/kg/s)\r
120 !-----------------------------------------------------------------------\r
121 !-----------------------------------------------------------------------\r
122 IMPLICIT NONE\r
123 \r
124 !.......Arguments\r
125 ! DECLARATIONS - INTEGER\r
126 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &\r
127 ims,ime, jms,jme, kms,kme, &\r
128 its,ite, jts,jte, kts,kte, XTIME\r
129 \r
130 ! DECLARATIONS - REAL\r
131 REAL, INTENT(IN) :: DTPBL, EP1, &\r
132 G, ROVCP, RD, CPD\r
133 \r
134 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &\r
135 INTENT(IN) :: U3D, V3D, &\r
136 PP3D, DZ8W, T3D, &\r
137 QV3D, QC3D, QI3D, &\r
138 RR3D, TH3D\r
139 \r
140 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: PSIM, GZ1OZ0, &\r
141 HFX, QFX, TSK, &\r
142 PSFC, WSPD, MUT\r
143 \r
144 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH, REGIME, & \r
145 UST, RMOL\r
146 \r
147 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &\r
148 INTENT(INOUT) :: RUBLTEN, RVBLTEN, &\r
149 RTHBLTEN, RQVBLTEN, &\r
150 RQCBLTEN, RQIBLTEN\r
151 \r
152 real, dimension( ims:ime, kms:kme, jms:jme ), &\r
153 intent(inout) :: exch_h,exch_m\r
154 \r
155 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: KPBL2D\r
156 \r
157 #if (WRF_CHEM == 1)\r
158 !... Chem\r
159 INTEGER, INTENT(IN ) :: nchem, kdvel, ndvel, num_vert_mix\r
160 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, nchem ), INTENT(INOUT) :: CHEM3D\r
161 REAL, DIMENSION( ims:ime, kdvel, jms:jme, ndvel ), INTENT(IN) :: VD3D\r
162 #endif\r
163 !... Local Variables\r
164 \r
165 !... Integer\r
166 INTEGER :: I, J, K, L\r
167 !... Real\r
168 REAL RDT\r
169 REAL, PARAMETER :: KARMAN = 0.4\r
170 \r
171 #if (WRF_CHEM == 1)\r
172 !... Chem\r
173 REAL, DIMENSION( ims:ime, kms:kme, nchem ) :: CHEM2D\r
174 REAL, DIMENSION( ims:ime, kdvel, ndvel ) :: VD2D\r
175 #endif\r
176 \r
177 RDT = 1.0 / DTPBL\r
178 \r
179 !==================================================\r
180 \r
181 DO j = jts,jte \r
182 #if (WRF_CHEM == 1)\r
183 DO L = 1, nchem\r
184 DO K = kms,kme\r
185 DO I = ims, ime\r
186 CHEM2D(i,k,l) = chem3d(i,k,j,l)\r
187 ENDDO\r
188 ENDDO\r
189 ENDDO\r
190 \r
191 DO L = 1, ndvel\r
192 DO K = 1, kdvel\r
193 DO I = ims, ime\r
194 VD2D(i,k,l) = VD3D(i,k,j,l)\r
195 ENDDO\r
196 ENDDO\r
197 ENDDO\r
198 #endif\r
199 CALL ACM2D(j=j,xtime=xtime,dtpbl=dtpbl &\r
200 ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j) &\r
201 ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j) &\r
202 ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j) &\r
203 ,qis=qi3d(ims,kms,j) &\r
204 #if (WRF_CHEM == 1)\r
205 ,chem=chem2d &\r
206 ,vd=vd2d &\r
207 ,nchem=nchem,kdvel=kdvel,ndvel=ndvel &\r
208 ,num_vert_mix=num_vert_mix &\r
209 #endif\r
210 ,dzf=DZ8W(ims,kms,j) &\r
211 ,densx=RR3D(ims,kms,j) &\r
212 ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &\r
213 ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j) &\r
214 ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &\r
215 ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt &\r
216 ,psfcpa=psfc(ims,j),ust=ust(ims,j) &\r
217 ,pbl=pblh(ims,j) &\r
218 ,exch_hx=exch_h(ims,kms,j) &\r
219 ,exch_mx=exch_m(ims,kms,j) &\r
220 ,regime=regime(ims,j),psim=psim(ims,j) &\r
221 ,hfx=hfx(ims,j),qfx=qfx(ims,j) &\r
222 ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &\r
223 ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j) &\r
224 ,mut=mut(ims,j), rmol=rmol(ims,j) &\r
225 ,ep1=ep1,karman=karman &\r
226 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &\r
227 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &\r
228 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )\r
229 #if (WRF_CHEM == 1)\r
230 DO L = 1, nchem\r
231 DO I = ims, ime\r
232 chem3d(i,kms:kme,j,l) = CHEM2D(i,kms:kme,l)\r
233 ENDDO\r
234 ENDDO\r
235 #endif\r
236 ENDDO\r
237 \r
238 END SUBROUTINE ACMPBL\r
239 !-----------------------------------------------------------------------\r
240 !-----------------------------------------------------------------------\r
241 \r
242 \r
243 !-----------------------------------------------------------------------\r
244 !-----------------------------------------------------------------------\r
245 SUBROUTINE ACM2D(j,xtime,dtpbl &\r
246 ,us,vs,theta,tt,qvs,qcs,qis &\r
247 #if (WRF_CHEM == 1)\r
248 ,chem, vd, nchem, kdvel, ndvel &\r
249 ,num_vert_mix &\r
250 #endif\r
251 ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp &\r
252 ,cpd,g,rovcp,rd,rdt,psfcpa,ust &\r
253 ,pbl,exch_hx,exch_mx,regime,psim &\r
254 ,hfx,qfx,tg,gz1oz0,wspd ,klpbl &\r
255 ,mut, rmol &\r
256 ,ep1,karman &\r
257 ,ids,ide, jds,jde, kds,kde &\r
258 ,ims,ime, jms,jme, kms,kme &\r
259 ,its,ite, jts,jte, kts,kte )\r
260 \r
261 !-----------------------------------------------------------------------\r
262 !-----------------------------------------------------------------------\r
263 IMPLICIT NONE\r
264 \r
265 !.......Arguments\r
266 \r
267 !... Real\r
268 REAL , INTENT(IN) :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT\r
269 REAL , DIMENSION( ims:ime ), INTENT(INOUT) :: PBL, UST\r
270 \r
271 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, TT, &\r
272 QVS, QCS, QIS, DENSX\r
273 REAL, DIMENSION( ims:ime, kms:kme ), intent(in) :: DZF\r
274 REAL, DIMENSION( ims:ime, kms:kme ), intent(inout) :: utnp, &\r
275 vtnp, &\r
276 ttnp, &\r
277 qvtnp, &\r
278 qctnp, &\r
279 qitnp\r
280 real, dimension( ims:ime ), intent(in ) :: psfcpa\r
281 real, dimension( ims:ime ), intent(in ) :: tg\r
282 real, dimension( ims:ime ), intent(inout) :: regime, rmol\r
283 real, dimension( ims:ime ), intent(in) :: wspd, psim, gz1oz0\r
284 real, dimension( ims:ime ), intent(in) :: hfx, qfx\r
285 real, dimension( ims:ime ), intent(in) :: mut\r
286 real, dimension( ims:ime, kms:kme ), &\r
287 intent(inout) :: exch_hx,exch_mx\r
288 !... Integer\r
289 INTEGER, DIMENSION( ims:ime ), INTENT(OUT):: KLPBL\r
290 INTEGER, INTENT(IN) :: XTIME\r
291 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &\r
292 ims,ime, jms,jme, kms,kme, &\r
293 its,ite, jts,jte, kts,kte, j\r
294 #if (WRF_CHEM == 1)\r
295 !....Chem\r
296 INTEGER, INTENT(IN) :: NCHEM, KDVEL,NDVEL,NUM_VERT_MIX\r
297 REAL , DIMENSION( ims:ime, kms:kme, NCHEM ), INTENT(INOUT) :: CHEM\r
298 REAL , DIMENSION( ims:ime, KDVEL, NDVEL ), INTENT(IN) :: VD\r
299 #endif\r
300 !--------------------------------------------------------------------\r
301 !--Local \r
302 INTEGER I, K \r
303 INTEGER :: KPBLHT\r
304 INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV, KSRC\r
305 \r
306 !... Real\r
307 REAL :: TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ\r
308 REAL :: ZH1,UH1,VH1 ! NEW FOR V3.7\r
309 REAL :: psix, THV1\r
310 REAL, DIMENSION( its:ite ) :: FINT, PSTAR, CPAIR\r
311 REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB, &\r
312 EDDYZ, EDDYZM, UX, VX, THETAX, &\r
313 QVX, QCX, QIX, ZA, DZHI, DZFI, DZH\r
314 REAL, DIMENSION( its:ite, 0:kte ) :: ZF\r
315 REAL, DIMENSION( its:ite) :: WST, TST, QST, USTM, TSTV\r
316 REAL, DIMENSION( its:ite ) :: MOL\r
317 REAL :: FINTT, ZMIX, UMIX, VMIX\r
318 REAL :: TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH\r
319 REAL :: A,TST12,RL,ZFUNC,DENSF\r
320 \r
321 !... Integer\r
322 INTEGER :: KL, jtf, ktf, itf, KMIX\r
323 character*512 :: message\r
324 !-----initialize vertical tendencies and\r
325 \r
326 DO i = its,ite\r
327 DO k = kts,kte\r
328 utnp(i,k) = 0.\r
329 vtnp(i,k) = 0.\r
330 ttnp(i,k) = 0.\r
331 ENDDO\r
332 ENDDO\r
333 \r
334 DO k = kts,kte\r
335 DO i = its,ite\r
336 qvtnp(i,k) = 0.\r
337 ENDDO\r
338 ENDDO\r
339 \r
340 DO k = kts,kte\r
341 DO i = its,ite\r
342 qctnp(i,k) = 0.\r
343 qitnp(i,k) = 0.\r
344 ENDDO\r
345 ENDDO\r
346 \r
347 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
348 !!! Compute Micromet Scaling variables, not availiable in WRF for ACM\r
349 DO I = its,ite\r
350 CPAIR(I) = CPD * (1.0 + 0.84 * QVS(I,1)) ! J/(K KG)\r
351 TMPFX = HFX(I) / (cpair(i) * DENSX(I,1))\r
352 TMPVTCON = 1.0 + EP1 * QVS(I,1) ! COnversion factor for virtual temperature\r
353 WS1 = SQRT(US(I,1)**2 + VS(I,1)**2) ! Level 1 wind speed\r
354 TST(I) = -TMPFX / UST(I)\r
355 QST(I) = -QFX(I) / (UST(I)*DENSX(I,1))\r
356 USTM(I) = UST(I) * WS1 / wspd(i)\r
357 THV1 = TMPVTCON * THETA(I,1) \r
358 TSTV(I) = TST(I)*TMPVTCON + THV1*EP1*QST(I)\r
359 IF(ABS(TSTV(I)).LT.1.0E-6) THEN\r
360 TSTV(I) = SIGN(1.0E-6,TSTV(I))\r
361 ENDIF\r
362 MOL(I) = THV1 * UST(i)**2/(KARMAN*G*TSTV(I))\r
363 RMOL(I) = 1.0/MOL(I)\r
364 WST(I) = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333 \r
365 PSTAR(I) = MUT(I)/1000. ! P* in cb \r
366 ENDDO\r
367 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
368 \r
369 !... Compute PBL height\r
370 \r
371 !... compute the height of full- and half-sigma level above ground level\r
372 DO I = its,ite\r
373 ZF(I,0) = 0.0\r
374 KLPBL(I) = 1\r
375 ENDDO\r
376 \r
377 ! write(0,*) 'HHHHH2'\r
378 DO K = kts, kte\r
379 DO I = its,ite\r
380 ZF(I,K) = DZF(I,K) + ZF(I,K-1)\r
381 ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))\r
382 DZH(I,K) = ZF(I,K) - ZF(I,K-1)\r
383 DZHI(I,K)= 1./DZH(I,K)\r
384 ENDDO\r
385 ENDDO\r
386 \r
387 DO K = kts, kte-1\r
388 DO I = its,ite\r
389 DZFI(I,K) = 1./(ZA(I,K+1)-ZA(I,K))\r
390 ENDDO\r
391 ENDDO\r
392 DO I = its,ite\r
393 DZFI(I,KTE) = DZFI(I,KTE-1)\r
394 ENDDO\r
395 \r
396 DO K = kts, kte\r
397 DO I = its,ite\r
398 TVCON = 1.0 + EP1 * QVS(I,K)\r
399 THETAV(I,K) = THETA(I,K) * TVCON\r
400 ENDDO\r
401 ENDDO\r
402 \r
403 \r
404 !... COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990 \r
405 DO 100 I = its,ite\r
406 DO K = 1,kte\r
407 RIB(I,K) = 0.0\r
408 KSRC(I) = K\r
409 IF (ZF(I,K).gt.30.0) GO TO 69\r
410 ENDDO\r
411 69 CONTINUE \r
412 \r
413 TH1 = 0.0\r
414 ZH1 = 0.0\r
415 UH1 = 0.0\r
416 VH1 = 0.0\r
417 DO K = 1,KSRC(I)\r
418 TH1 = TH1 + THETAV(I,K) \r
419 ZH1 = ZH1 + ZA(I,K)\r
420 UH1 = UH1 + US(I,K)\r
421 VH1 = VH1 + VS(I,K)\r
422 ENDDO \r
423 TH1 = TH1/KSRC(I)\r
424 ZH1 = ZH1/KSRC(I)\r
425 UH1 = UH1/KSRC(I)\r
426 VH1 = VH1/KSRC(I)\r
427 IF(MOL(I).LT.0.0 .AND. XTIME.GT.1) then\r
428 WSS = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333\r
429 TCONV = -8.5 * UST(I) * TSTV(I) / WSS\r
430 TH1 = TH1 + TCONV\r
431 ENDIF\r
432 \r
433 99 KMIX = KSRC(I)\r
434 DO K = KSRC(I),kte\r
435 DTMP = THETAV(I,K) - TH1\r
436 IF (DTMP.LT.0.0) KMIX = K\r
437 ENDDO\r
438 IF(KMIX.GT.KSRC(I)) THEN\r
439 FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1) &\r
440 - THETAV(I,KMIX))\r
441 ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)\r
442 UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)\r
443 VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)\r
444 ELSE\r
445 ZMIX = ZH1\r
446 UMIX = UH1\r
447 VMIX = VH1\r
448 ENDIF\r
449 DO K = KMIX,kte\r
450 DTMP = THETAV(I,K) - TH1\r
451 TOG = 0.5 * (THETAV(I,K) + TH1) / G\r
452 WSSQ = (US(I,K)-UMIX)**2 &\r
453 + (VS(I,K)-VMIX)**2\r
454 IF (KMIX .EQ. KSRC(I)) THEN\r
455 WSSQ = US(I,K)**2 + VS(I,K)**2\r
456 WSSQ = WSSQ + 100.*UST(I)*UST(I) \r
457 ENDIF\r
458 WSSQ = MAX( WSSQ, 0.1 )\r
459 RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)\r
460 IF (RIB(I,K) .GE. RIC) GO TO 201\r
461 ENDDO\r
462 \r
463 \r
464 write (message,*)' RIBX never exceeds RIC, RIB(i,kte) = ',rib(i,5), &\r
465 ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i), &\r
466 ' TCONV = ',TCONV,' WST = ',WST(I), &\r
467 ' KMIX = ',kmix,' UST = ',UST(I), &\r
468 ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1), &\r
469 ' I,J=',I,J\r
470 CALL wrf_error_fatal ( message )\r
471 \r
472 201 CONTINUE\r
473 \r
474 KPBLH(I) = K\r
475 \r
476 100 CONTINUE\r
477 \r
478 DO I = its,ite\r
479 IF (KPBLH(I) .GT. KSRC(I)) THEN\r
480 !---------INTERPOLATE BETWEEN LEVELS -- jp 7/93\r
481 FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) - &\r
482 RIB(I,KPBLH(I)-1))\r
483 IF (FINT(I) .GT. 0.5) THEN\r
484 KPBLHT = KPBLH(I)\r
485 FINT(I) = FINT(I) - 0.5\r
486 ELSE\r
487 KPBLHT = KPBLH(I) - 1\r
488 FINT(I) = FINT(I) + 0.5\r
489 ENDIF\r
490 PBL(I) = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) + &\r
491 ZF(I,KPBLHT-1)\r
492 KLPBL(I) = KPBLHT\r
493 ELSE\r
494 KLPBL(I) = KSRC(I)\r
495 PBL(I) = ZA(I,KSRC(I)) \r
496 ENDIF\r
497 \r
498 ENDDO\r
499 \r
500 DO I = its,ite \r
501 NOCONV(I) = 0\r
502 \r
503 ! Check for CBL and identify conv. vs. non-conv cells\r
504 IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3 &\r
505 .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN\r
506 NOCONV(I) = 1\r
507 REGIME(I) = 4.0 ! FREE CONVECTIVE - ACM\r
508 ENDIF\r
509 ENDDO\r
510 \r
511 !... Calculate Kz\r
512 CALL EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &\r
513 US, VS, TT, THETAV, DENSX, PSTAR, &\r
514 QVS, QCS, QIS, G, RD, CPAIR, &\r
515 EDDYZ, EDDYZM, its,ite, kts,kte,ims,ime, kms,kme)\r
516 \r
517 CALL ACM(DTPBL, PSTAR, NOCONV, ZF, DZH, DZHI, J, &\r
518 KLPBL, PBL, DZFI, MOL, UST, &\r
519 TST, QST, USTM, EDDYZ, DENSX, &\r
520 THETA, QVS, QCS, QIS, &\r
521 THETAX, QVX, QCX, QIX, &\r
522 #if (WRF_CHEM == 1)\r
523 CHEM, VD, NCHEM, KDVEL, NDVEL,NUM_VERT_MIX, &\r
524 #endif\r
525 ids,ide, jds,jde, kds,kde, &\r
526 ims,ime, jms,jme, kms,kme, &\r
527 its,ite, jts,jte, kts,kte)\r
528 CALL ACMM(DTPBL, PSTAR, NOCONV, ZF, DZH, DZHI, J, &\r
529 KLPBL, PBL, DZFI, MOL, UST, &\r
530 TST, QST, USTM, EDDYZM, DENSX, &\r
531 US, VS, &\r
532 UX, VX, &\r
533 ids,ide, jds,jde, kds,kde, &\r
534 ims,ime, jms,jme, kms,kme, &\r
535 its,ite, jts,jte, kts,kte)\r
536 \r
537 !.. Load exch_h for use in CCN activation\r
538 \r
539 DO K = kts, kte-1\r
540 DO I = its, ite\r
541 exch_hx(I,K) = EDDYZ(I,K)\r
542 exch_mx(I,K) = EDDYZM(I,K)\r
543 ENDDO\r
544 ENDDO\r
545 \r
546 !... Calculate tendency due to PBL parameterization\r
547 \r
548 DO K = kts, kte\r
549 DO I = its, ite\r
550 UTNP(I,K) = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT\r
551 VTNP(I,K) = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT\r
552 TTNP(I,K) = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT\r
553 QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT\r
554 QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT\r
555 QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT\r
556 ENDDO\r
557 ENDDO\r
558 \r
559 END SUBROUTINE ACM2D\r
560 !-----------------------------------------------------------------------\r
561 !-----------------------------------------------------------------------\r
562 \r
563 !-----------------------------------------------------------------------\r
564 !-----------------------------------------------------------------------\r
565 SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &\r
566 RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &\r
567 restart, allowed_to_read , &\r
568 ids, ide, jds, jde, kds, kde, &\r
569 ims, ime, jms, jme, kms, kme, &\r
570 its, ite, jts, jte, kts, kte )\r
571 !-----------------------------------------------------------------------\r
572 !\r
573 ! This subroutine is for preparing ACM PBL variables. \r
574 ! Called from module_physics_init.F\r
575 !\r
576 ! REVISION HISTORY:\r
577 ! AX 3/2005 - Originally developed\r
578 !-----------------------------------------------------------------------\r
579 ! ARGUMENT LIST:\r
580 !\r
581 !-----------------------------------------------------------------------\r
582 !-----------------------------------------------------------------------\r
583 IMPLICIT NONE\r
584 !\r
585 LOGICAL , INTENT(IN) :: restart , allowed_to_read\r
586 \r
587 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &\r
588 ims, ime, jms, jme, kms, kme, &\r
589 its, ite, jts, jte, kts, kte\r
590 \r
591 INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR\r
592 \r
593 ! REAL , DIMENSION( kms:kme ), INTENT(IN) :: SHALF\r
594 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &\r
595 RUBLTEN, &\r
596 RVBLTEN, &\r
597 RTHBLTEN, &\r
598 RQVBLTEN, &\r
599 RQCBLTEN, & \r
600 RQIBLTEN\r
601 \r
602 !... Local Variables\r
603 INTEGER :: i, j, k, itf, jtf, ktf\r
604 \r
605 !\r
606 jtf=min0(jte,jde-1)\r
607 ktf=min0(kte,kde-1)\r
608 itf=min0(ite,ide-1)\r
609 \r
610 IF(.not.restart)THEN\r
611 DO j=jts,jtf\r
612 DO k=kts,ktf\r
613 DO i=its,itf\r
614 RUBLTEN(i,k,j)=0.\r
615 RVBLTEN(i,k,j)=0.\r
616 RTHBLTEN(i,k,j)=0.\r
617 RQVBLTEN(i,k,j)=0.\r
618 RQCBLTEN(i,k,j)=0.\r
619 ENDDO\r
620 ENDDO\r
621 ENDDO\r
622 ENDIF\r
623 \r
624 IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN\r
625 DO j=jts,jtf\r
626 DO k=kts,ktf\r
627 DO i=its,itf\r
628 RQIBLTEN(i,k,j)=0.\r
629 ENDDO\r
630 ENDDO\r
631 ENDDO\r
632 ENDIF\r
633 \r
634 \r
635 END SUBROUTINE acminit\r
636 !-----------------------------------------------------------------------\r
637 !-----------------------------------------------------------------------\r
638 \r
639 !-----------------------------------------------------------------------\r
640 !------------------------------------------------------------------- \r
641 SUBROUTINE EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &\r
642 US, VS, TT, THETAV, DENSX, PSTAR, &\r
643 QVS, QCS, QIS, G, RD, CPAIR, &\r
644 EDDYZ, EDDYZM, its,ite, kts,kte,ims,ime,kms,kme )\r
645 \r
646 \r
647 !**********************************************************************\r
648 ! Two methods for computing Kz:\r
649 ! 1. Boundary scaling similar to Holtslag and Boville (1993)\r
650 ! 2. Local Kz computed as function of local Richardson # and vertical \r
651 ! wind shear, similar to LIU & CARROLL (1996)\r
652 !\r
653 !**********************************************************************\r
654 !\r
655 !-- DTPBL time step of the minor loop for the land-surface/pbl model\r
656 !-- ZF height of full level\r
657 !-- ZA height of half level\r
658 !-- MOL Monin-Obukhov length in 1D form\r
659 !-- PBL PBL height in 1D form\r
660 !-- UST friction velocity U* in 1D form (m/s)\r
661 !-- US U wind \r
662 !-- VS V wind\r
663 !-- TT temperature\r
664 !-- THETAV potential virtual temperature\r
665 !-- DENSX dry air density (kg/m^3)\r
666 !-- PSTAR P*=Psfc-Ptop\r
667 !-- QVS water vapor mixing ratio (Kg/Kg)\r
668 !-- QCS cloud mixing ratio (Kg/Kg)\r
669 !-- QIS ice mixing ratio (Kg/Kg)\r
670 !-- G gravity\r
671 !-- RD gas constant for dry air (j/kg/k)\r
672 !-- CPAIR specific heat of moist air (M^2 S^-2 K^-1)\r
673 !-- EDDYZ eddy diffusivity for heat KZ\r
674 !-- EDDYZM eddy diffusivity for momentum KM\r
675 !-----------------------------------------------------------------------\r
676 !-----------------------------------------------------------------------\r
677 \r
678 IMPLICIT NONE\r
679 \r
680 !.......Arguments\r
681 \r
682 !... Integer\r
683 INTEGER, INTENT(IN) :: its,ite, kts,kte,ims,ime,kms,kme\r
684 !... Real\r
685 REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST\r
686 REAL , INTENT(IN) :: DTPBL, G, RD\r
687 REAL , DIMENSION( its:ite ), INTENT(IN) :: MOL, PSTAR, CPAIR\r
688 \r
689 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, TT, &\r
690 QVS, QCS, QIS, DENSX\r
691 REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV\r
692 REAL, DIMENSION( its:ite, 0:kte ) , INTENT(IN) :: ZF\r
693 \r
694 REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ,EDDYZM\r
695 \r
696 !.......Local variables\r
697 \r
698 !... Integer\r
699 INTEGER :: ILX, KL, KLM, K, I\r
700 \r
701 !... Real\r
702 REAL :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ\r
703 REAL :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF, KZO\r
704 REAL :: FH, FM\r
705 REAL :: WM, EDYZM, PHIM\r
706 !... Parameters\r
707 REAL, PARAMETER :: RV = 461.5\r
708 REAL, PARAMETER :: RC = 0.25\r
709 REAL, PARAMETER :: RLAM = 80.0\r
710 REAL, PARAMETER :: GAMH = 16.0 !Dyer74 !15.0 ! Holtslag and Boville (1993)\r
711 REAL, PARAMETER :: GAMM = 16.0 !Dyer74\r
712 REAL, PARAMETER :: BETAH = 5.0 ! Holtslag and Boville (1993) BETAM = BETAH\r
713 REAL, PARAMETER :: KARMAN = 0.4\r
714 REAL, PARAMETER :: P = 2.0 ! ZFUNC exponent\r
715 REAL, PARAMETER :: EDYZ0 = 0.01 ! New Min Kz\r
716 REAL, PARAMETER :: PR = 0.8 ! Prandtl #\r
717 INTEGER, PARAMETER :: imvdif = 1\r
718 !\r
719 ILX = ite \r
720 KL = kte\r
721 KLM = kte - 1\r
722 \r
723 DO K = kts,KLM\r
724 DO I = its,ILX\r
725 EDYZ = 0.0\r
726 ZOVL = 0.0\r
727 DZF = ZA(I,K+1) - ZA(I,K)\r
728 KZO = EDYZ0\r
729 !--------------------------------------------------------------------------\r
730 IF (ZF(I,K) .LT. PBL(I)) THEN\r
731 ZOVL = ZF(I,K) / MOL(I)\r
732 IF (ZOVL .LT. 0.0) THEN\r
733 IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN\r
734 PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)\r
735 PHIM = (1.0 - GAMM * ZOVL)**(-0.25)\r
736 WT = UST(I) / PHIH\r
737 WM = UST(I) / PHIM\r
738 ELSE\r
739 ZSOL = 0.1 * PBL(I) / MOL(I)\r
740 PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)\r
741 PHIM = (1.0 - GAMM * ZSOL)**(-0.25)\r
742 WT = UST(I) / PHIH\r
743 WM = UST(I) / PHIM\r
744 ENDIF\r
745 ELSE IF (ZOVL .LT. 1.0) THEN\r
746 PHIH = 1.0 + BETAH * ZOVL\r
747 WT = UST(I) / PHIH\r
748 WM = WT\r
749 ELSE\r
750 PHIH = BETAH + ZOVL\r
751 WT = UST(I) / PHIH\r
752 WM = WT\r
753 ENDIF\r
754 ZFUNC = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** P\r
755 EDYZ = KARMAN * WT * ZFUNC\r
756 EDYZM = KARMAN * WM * ZFUNC\r
757 ENDIF\r
758 !--------------------------------------------------------------------------!--------------------------------------------------------------------------\r
759 SS = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2) &\r
760 / (DZF * DZF) + 1.0E-9\r
761 GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))\r
762 RI = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)\r
763 !--------------------------------------------------------------------------\r
764 ! Adjustment to vert diff in Moist air\r
765 IF(imvdif.eq.1)then\r
766 IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+ &\r
767 QIS(I,K+1)) .GT. 0.01E-3) THEN\r
768 QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))\r
769 TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))\r
770 XLV = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6\r
771 ALPH = XLV * QMEAN / RD / TMEAN\r
772 CHI = XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN\r
773 RI = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) * &\r
774 ((CHI - ALPH) / (1.0 + CHI)))\r
775 ENDIF\r
776 ENDIF\r
777 !--------------------------------------------------------------------------\r
778 \r
779 ZK = 0.4 * ZF(I,K)\r
780 SQL = (ZK * RLAM / (RLAM + ZK)) ** 2\r
781 \r
782 IF (RI .GE. 0.0) THEN\r
783 FH=1./(1.+10.*RI+50.*RI**2+5000.*RI**4)+0.0012 !pleim5\r
784 FM= PR*FH + 0.00104\r
785 \r
786 EDDYZ(I,K) = KZO + SQRT(SS) * FH * SQL\r
787 EDDYZM(I,K) = KZO + SQRT(SS) * FM * SQL\r
788 ELSE\r
789 EDDYZ(I,K) = KZO + SQRT(SS * (1.0 - 25.0 * RI)) * SQL\r
790 EDDYZM(I,K) = EDDYZ(I,K) * PR\r
791 ENDIF\r
792 \r
793 IF(EDYZ.GT.EDDYZ(I,K)) THEN\r
794 EDDYZ(I,K) = EDYZ\r
795 EDDYZM(I,K) = MIN(EDYZM,EDYZ*0.8) !PR\r
796 ENDIF\r
797 \r
798 EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))\r
799 EDDYZ(I,K) = MAX(KZO,EDDYZ(I,K))\r
800 EDDYZM(I,K) = MIN(1000.0,EDDYZM(I,K))\r
801 EDDYZM(I,K) = MAX(KZO,EDDYZM(I,K))\r
802 \r
803 ENDDO ! for I loop\r
804 ENDDO ! for k loop\r
805 !\r
806 DO I = its,ILX\r
807 EDDYZ(I,KL) = 0.0\r
808 EDDYZM(I,KL) = 0.0\r
809 ENDDO\r
810 \r
811 END SUBROUTINE EDDYX\r
812 !-----------------------------------------------------------------------\r
813 !-----------------------------------------------------------------------\r
814 \r
815 !-----------------------------------------------------------------------\r
816 !------------------------------------------------------------------- \r
817 SUBROUTINE ACM (DTPBL, PSTAR, NOCONV, ZF, DZH, DZHI, JX, &\r
818 KLPBL, PBL, DZFI, MOL, UST, &\r
819 TST, QST, USTM, EDDYZ, DENSX, &\r
820 THETA, QVS, QCS, QIS, &\r
821 THETAX, QVX, QCX, QIX, &\r
822 #if (WRF_CHEM == 1)\r
823 CHEM, VD, NCHEM, KDVEL, NDVEL, &\r
824 NUM_VERT_MIX, &\r
825 #endif\r
826 ids,ide, jds,jde, kds,kde, &\r
827 ims,ime, jms,jme, kms,kme, &\r
828 its,ite, jts,jte, kts,kte)\r
829 !**********************************************************************\r
830 ! PBL model called the Asymmetric Convective Model, Version 2 (ACM2) \r
831 ! -- See top of module for summary and references\r
832 !\r
833 !---- REVISION HISTORY:\r
834 ! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM\r
835 ! JP and RG 8/2006 - updates\r
836 ! JP 3/2013 - Chem additions\r
837 !\r
838 !**********************************************************************\r
839 ! ARGUMENTS:\r
840 !-- DTPBL PBL time step\r
841 !-- PSTAR Psurf - Ptop in cb\r
842 !-- NOCONV If free convection =0, no; =1, yes\r
843 !-- ZF Height for full layer\r
844 !-- DZH Layer thickness --> ZF(K) - ZF(K-1)\r
845 !-- DZHI Inverse of layer thickness\r
846 !-- JX N-S index\r
847 !-- KLPBL PBL level at K index\r
848 !-- PBL PBL height in m\r
849 !-- DZFI Inverse layer thickness --> 1/(Z(K+1)-Z(K))\r
850 !-- MOL Monin-Obukhov length in 1D form\r
851 !-- UST U* in 1D form\r
852 !-- TST Theta* in 1D form\r
853 !-- QST Q* in 1D form\r
854 !-- USTM U* for computation of momemtum flux \r
855 !-- EDDYZ eddy diffusivity KZ\r
856 !-- DENSX dry air density (kg/m^3)\r
857 !-- US U wind \r
858 !-- VS V wind\r
859 !-- THETA potential temperature\r
860 !-- QVS water vapor mixing ratio (Kg/Kg)\r
861 !-- QCS cloud mixing ratio (Kg/Kg)\r
862 !-- QIS ice mixing ratio (Kg/Kg)\r
863 !-- UX new U wind \r
864 !-- VX new V wind\r
865 !-- THETAX new potential temperature\r
866 !-- QVX new water vapor mixing ratio (Kg/Kg)\r
867 !-- QCX new cloud mixing ratio (Kg/Kg)\r
868 !-- QIX new ice mixing ratio (Kg/Kg)\r
869 !-- CHEM Chemical species mixing ratios (ppm) WRFChem \r
870 !-- VD Dry deposition velocity (m/s) WRFChem\r
871 !-- NCHEM Number of chemical species WRFChem\r
872 !-----------------------------------------------------------------------\r
873 !-----------------------------------------------------------------------\r
874 \r
875 IMPLICIT NONE\r
876 \r
877 !.......Arguments\r
878 \r
879 !... Integer\r
880 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &\r
881 ims,ime, jms,jme, kms,kme, &\r
882 its,ite, jts,jte, kts,kte, JX\r
883 INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV\r
884 INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL\r
885 \r
886 !... Real\r
887 REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST\r
888 REAL , INTENT(IN) :: DTPBL\r
889 REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, &\r
890 MOL, TST, &\r
891 QST, USTM\r
892 REAL , DIMENSION( its:ite, kts:kte ), INTENT(IN) :: DZHI, DZH, DZFI\r
893 REAL , DIMENSION( its:ite, 0:kte ), INTENT(IN) :: ZF\r
894 REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZ\r
895 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: THETA, &\r
896 QVS, QCS, QIS, DENSX\r
897 REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: THETAX, &\r
898 QVX, QCX, QIX\r
899 #if (WRF_CHEM == 1)\r
900 !......Chem\r
901 INTEGER, INTENT(IN) :: NCHEM, KDVEL, NDVEL, NUM_VERT_MIX\r
902 REAL , DIMENSION( ims:ime, kms:kme, nchem ), INTENT(INOUT) :: CHEM \r
903 REAL , DIMENSION( ims:ime, KDVEL, NDVEL ), INTENT(IN) :: VD\r
904 #endif\r
905 !.......Local variables\r
906 \r
907 !... Parameters\r
908 INTEGER, PARAMETER :: NSP = 4 \r
909 !\r
910 REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP\r
911 REAL, PARAMETER :: KARMAN = 0.4\r
912 \r
913 !... Integer\r
914 INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L,LL\r
915 INTEGER :: KCBLMX\r
916 INTEGER, DIMENSION( its:ite ) :: KCBL\r
917 \r
918 !... Real\r
919 REAL :: MBMAX, HOVL, MEDDY, MBAR\r
920 REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1\r
921 REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, DTLIM\r
922 REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN\r
923 REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS\r
924 REAL DELC\r
925 REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI\r
926 REAL, ALLOCATABLE, DIMENSION( : , : ) :: DI, UI \r
927 REAL, ALLOCATABLE, DIMENSION( : , : ) :: FS\r
928 REAL, ALLOCATABLE, DIMENSION( : , : , : ) :: VCI\r
929 \r
930 CHARACTER*80 :: message\r
931 \r
932 !\r
933 !--Start Exicutable ----\r
934 \r
935 ILX = ite\r
936 KL = kte\r
937 KLM = kte - 1\r
938 NSPX = NSP\r
939 #if (WRF_CHEM == 1)\r
940 NSPX = NSPX + NUM_VERT_MIX \r
941 #endif\r
942 \r
943 KCBLMX = 0\r
944 MBMAX = 0.0\r
945 !...Allocate species variables\r
946 ALLOCATE (DI( 1:NSPX,kts:kte )) \r
947 ALLOCATE (UI( 1:NSPX,kts:kte )) \r
948 ALLOCATE (FS( 1:NSPX, its:ite )) \r
949 ALLOCATE (VCI( 1:NSPX,its:ite,kts:kte ))\r
950 \r
951 !---COMPUTE ACM MIXING RATE\r
952 DO I = its, ILX\r
953 DTLIM(I) = DTPBL\r
954 PSTARI(I) = 1.0 / PSTAR(I)\r
955 KCBL(I) = 1\r
956 FSACM(I) = 0.0\r
957 \r
958 IF (NOCONV(I) .EQ. 1) THEN\r
959 KCBL(I) = KLPBL(I)\r
960 \r
961 !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE\r
962 !--New couple ACM & EDDY-------------------------------------------------------------\r
963 HOVL = -PBL(I) / MOL(I)\r
964 FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))\r
965 MEDDY = EDDYZ(I,1) * DZFI(i,1) / (PBL(I) - ZF(i,1))\r
966 MBAR = MEDDY * FSACM(I)\r
967 DO K = kts,KCBL(I)-1\r
968 EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))\r
969 ENDDO\r
970 \r
971 MBMAX = AMAX1(MBMAX,MBAR)\r
972 DO K = kts+1,KCBL(I)\r
973 MBARKS(K,I) = MBAR\r
974 MDWN(K,I) = MBAR * (PBL(I) - ZF(i,K-1)) * DZHI(i,K)\r
975 ENDDO\r
976 MBARKS(1,I) = MBAR\r
977 MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)\r
978 MDWN(KCBL(I)+1,I) = 0.0\r
979 ENDIF\r
980 ENDDO ! end of I loop\r
981 \r
982 DO K = kts,KLM\r
983 DO I = its,ILX\r
984 EKZ = EDDYZ(I,K) * DZFI(i,K) * DZHI(i,K)\r
985 DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))\r
986 ENDDO\r
987 ENDDO\r
988 \r
989 DO I = its,ILX \r
990 IF (NOCONV(I) .EQ. 1) THEN\r
991 KCBLMX = AMAX0(KLPBL(I),KCBLMX)\r
992 RZ = (ZF(i,KCBL(I)) - ZF(i,1)) * DZHI(i,1)\r
993 DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))\r
994 ENDIF\r
995 ENDDO\r
996 \r
997 DO K = kts,KL\r
998 DO I = its,ILX\r
999 VCI(1,I,K) = THETA(I,K)\r
1000 VCI(2,I,K) = QVS(I,K)\r
1001 ! -- Also mix cloud water and ice IF necessary\r
1002 ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN !!! Check other PBL models\r
1003 VCI(3,I,K) = QCS(I,K)\r
1004 VCI(4,I,K) = QIS(I,K)\r
1005 #if (WRF_CHEM == 1)\r
1006 DO L= NSP+1, NSPX\r
1007 VCI(L,I,K) = CHEM(I,K,L-NSP)\r
1008 ENDDO\r
1009 #endif\r
1010 ENDDO\r
1011 ENDDO\r
1012 \r
1013 DO I = its,ILX\r
1014 FS(1,I) = -UST(I) * TST(I)\r
1015 FS(2,I) = -UST(I) * QST(I)\r
1016 FS(3,I) = 0.0\r
1017 FS(4,I) = 0.0 ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0\r
1018 #if (WRF_CHEM == 1)\r
1019 DO L= NSP+1, NSPX\r
1020 FS(L,I) = -VD(I,1,L-NSP) * CHEM(I,1,L-NSP)\r
1021 ENDDO\r
1022 #endif\r
1023 ENDDO\r
1024 !\r
1025 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
1026 DO I = its,ILX \r
1027 \r
1028 NLP = INT(DTPBL / DTLIM(I) + 1.0)\r
1029 DTS = (DTPBL / NLP)\r
1030 DTRAT = DTS / DTPBL\r
1031 DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP \r
1032 \r
1033 !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES\r
1034 \r
1035 DO K = kts,kte\r
1036 AI(K) = 0.0\r
1037 BI(K) = 0.0\r
1038 CI(K) = 0.0\r
1039 EI(K) = 0.0\r
1040 ENDDO\r
1041 \r
1042 DO K = 2, KCBL(I)\r
1043 EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DZH(i,K) * DZHI(i,K-1)\r
1044 BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS\r
1045 AI(K) = -CRANKP * MBARKS(K,I) * DTS\r
1046 ENDDO\r
1047 \r
1048 EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DZHI(i,1) * DZFI(i,1) * DTS\r
1049 AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DZHI(i,2) * DZFI(i,1) * DTS\r
1050 \r
1051 DO K = KCBL(I)+1, KL\r
1052 BI(K) = 1.0\r
1053 ENDDO\r
1054 \r
1055 DO K = 2,KL\r
1056 XPLUS(K) = EDDYZ(I,K) * DZHI(i,K) * DZFI(i,K) * DTS\r
1057 XMINUS(K) = EDDYZ(I,K-1) * DZHI(i,K) * DZFI(i,K-1) * DTS\r
1058 CI(K) = - XMINUS(K) * CRANKP\r
1059 EI(K) = EI(K) - XPLUS(K) * CRANKP\r
1060 BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP\r
1061 ENDDO\r
1062 \r
1063 IF (NOCONV(I) .EQ. 1) THEN\r
1064 BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBL(I) - ZF(i,1)) * DTS &\r
1065 * DZHI(i,1) + EDDYZ(I,1) * CRANKP * DZHI(i,1) * DZFI(i,1) * DTS\r
1066 ELSE\r
1067 BI(1) = 1.0 + EDDYZ(I,1) * CRANKP * DZHI(i,1) * DZFI(i,1) * DTS\r
1068 ENDIF\r
1069 \r
1070 \r
1071 DO K = 1,KL\r
1072 DO L = 1,NSPX \r
1073 DI(L,K) = 0.0\r
1074 ENDDO\r
1075 ENDDO\r
1076 !\r
1077 !** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION\r
1078 DO K = 2,KCBL(I)\r
1079 DO L = 1,NSPX \r
1080 DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * &\r
1081 VCI(L,I,K) + DZH(i,K+1) * DZHI(i,K) * &\r
1082 MDWN(K+1,I) * VCI(L,I,K+1))\r
1083 DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC\r
1084 ENDDO\r
1085 ENDDO\r
1086 \r
1087 DO K = KCBL(I)+1, KL\r
1088 DO L = 1,NSPX \r
1089 DI(L,K) = VCI(L,I,K)\r
1090 ENDDO\r
1091 ENDDO\r
1092 \r
1093 DO K = 2,KL\r
1094 IF (K .EQ. KL) THEN\r
1095 DO L = 1,NSPX \r
1096 DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * &\r
1097 (VCI(L,I,K) - VCI(L,I,K-1))\r
1098 ENDDO\r
1099 ELSE\r
1100 DO L = 1,NSPX \r
1101 DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * &\r
1102 (VCI(L,I,K+1) - VCI(L,I,K)) - &\r
1103 (1.0 - CRANKP) * XMINUS(K) * &\r
1104 (VCI(L,I,K) - VCI(L,I,K-1))\r
1105 ENDDO\r
1106 ENDIF\r
1107 ENDDO\r
1108 \r
1109 IF (NOCONV(I) .EQ. 1) THEN\r
1110 DO L = 1,NSPX \r
1111 DI(L,1) = VCI(L,I,1) + (FS(L,I) - (1.0 - CRANKP) &\r
1112 * (MBARKS(1,I) * &\r
1113 (PBL(I) - ZF(i,1)) * VCI(L,I,1) - &\r
1114 MDWN(2,I) * VCI(L,I,2) * DZH(i,2))) * DZHI(i,1) * DTS\r
1115 ENDDO\r
1116 ELSE\r
1117 DO L = 1,NSPX \r
1118 DI(L,1) = VCI(L,I,1) + FS(L,I) * DZHI(i,1) * DTS\r
1119 ENDDO\r
1120 ENDIF\r
1121 DO L = 1,NSPX \r
1122 DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DZHI(i,1) &\r
1123 * DZFI(i,1) * DTS * (VCI(L,I,2) - VCI(L,I,1))\r
1124 ENDDO\r
1125 IF ( NOCONV(I) .EQ. 1 ) THEN\r
1126 CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)\r
1127 ELSE\r
1128 CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)\r
1129 END IF\r
1130 !\r
1131 !-- COMPUTE NEW THETAV AND Q\r
1132 DO K = 1,KL\r
1133 DO L = 1,NSPX \r
1134 VCI(L,I,K) = UI(L,K)\r
1135 ENDDO\r
1136 ENDDO\r
1137 \r
1138 ENDDO ! END I LOOP\r
1139 ENDDO ! END SUB TIME LOOP\r
1140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
1141 \r
1142 !\r
1143 DO K = kts,KL\r
1144 DO I = its,ILX\r
1145 THETAX(I,K) = VCI(1,I,K)\r
1146 QVX(I,K) = VCI(2,I,K)\r
1147 QCX(I,K) = VCI(3,I,K)\r
1148 QIX(I,K) = VCI(4,I,K)\r
1149 #if (WRF_CHEM == 1)\r
1150 DO LL= 7, NSPX\r
1151 CHEM(I,K,LL-NSP) = VCI(LL,I,K) \r
1152 ENDDO\r
1153 #endif\r
1154 ENDDO\r
1155 ENDDO\r
1156 \r
1157 DEALLOCATE (DI) \r
1158 DEALLOCATE (UI) \r
1159 DEALLOCATE (FS)\r
1160 DEALLOCATE (VCI)\r
1161 \r
1162 END SUBROUTINE ACM\r
1163 !-----------------------------------------------------------------------\r
1164 !-----------------------------------------------------------------------\r
1165 !------------------------------------------------------------------- \r
1166 SUBROUTINE ACMM (DTPBL, PSTAR, NOCONV, ZF, DZH, DZHI, JX, &\r
1167 KLPBL, PBL, DZFI, MOL, UST, &\r
1168 TST, QST, USTM, EDDYZM, DENSX, &\r
1169 US, VS, &\r
1170 UX, VX, &\r
1171 ids,ide, jds,jde, kds,kde, &\r
1172 ims,ime, jms,jme, kms,kme, &\r
1173 its,ite, jts,jte, kts,kte)\r
1174 !**********************************************************************\r
1175 ! PBL model called the Asymmetric Convective Model, Version 2 (ACM2) \r
1176 ! -- See top of module for summary and references\r
1177 !\r
1178 !---- REVISION HISTORY:\r
1179 ! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM\r
1180 ! JP and RG 8/2006 - updates\r
1181 !\r
1182 !**********************************************************************\r
1183 ! ARGUMENTS:\r
1184 !-- DTPBL PBL time step\r
1185 !-- PSTAR Psurf - Ptop in cb\r
1186 !-- NOCONV If free convection =0, no; =1, yes\r
1187 !-- ZF Height for full layer\r
1188 !-- DZH Layer thickness --> ZF(K) - ZF(K-1)\r
1189 !-- DZHI Inverse of layer thickness\r
1190 !-- JX N-S index\r
1191 !-- KLPBL PBL level at K index\r
1192 !-- PBL PBL height in m\r
1193 !-- DZFI Inverse layer thickness --> 1/(Z(K+1)-Z(K))\r
1194 !-- MOL Monin-Obukhov length in 1D form\r
1195 !-- UST U* in 1D form\r
1196 !-- TST Theta* in 1D form\r
1197 !-- QST Q* in 1D form\r
1198 !-- USTM U* for computation of momemtum flux \r
1199 !-- EDDYZM eddy diffusivity for momentum KM\r
1200 !-- DENSX dry air density (kg/m^3)\r
1201 !-- US U wind \r
1202 !-- VS V wind\r
1203 !-- THETA potential temperature\r
1204 !-- QVS water vapor mixing ratio (Kg/Kg)\r
1205 !-- QCS cloud mixing ratio (Kg/Kg)\r
1206 !-- QIS ice mixing ratio (Kg/Kg)\r
1207 !-- UX new U wind \r
1208 !-- VX new V wind\r
1209 !-- THETAX new potential temperature\r
1210 !-- QVX new water vapor mixing ratio (Kg/Kg)\r
1211 !-- QCX new cloud mixing ratio (Kg/Kg)\r
1212 !-- QIX new ice mixing ratio (Kg/Kg)\r
1213 !-----------------------------------------------------------------------\r
1214 !-----------------------------------------------------------------------\r
1215 \r
1216 IMPLICIT NONE\r
1217 \r
1218 !.......Arguments\r
1219 \r
1220 !... Integer\r
1221 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &\r
1222 ims,ime, jms,jme, kms,kme, &\r
1223 its,ite, jts,jte, kts,kte, JX\r
1224 INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV\r
1225 INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL\r
1226 \r
1227 !... Real\r
1228 REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST\r
1229 REAL , INTENT(IN) :: DTPBL\r
1230 REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, &\r
1231 MOL, TST, &\r
1232 QST, USTM\r
1233 REAL , DIMENSION( its:ite, kts:kte ), INTENT(IN) :: DZHI, DZH, DZFI\r
1234 REAL , DIMENSION( its:ite, 0:kte ), INTENT(IN) :: ZF\r
1235 REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZM\r
1236 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, &\r
1237 DENSX\r
1238 REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX\r
1239 !.......Local variables\r
1240 \r
1241 !... Parameters\r
1242 INTEGER, PARAMETER :: NSP = 2\r
1243 !\r
1244 REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP\r
1245 REAL, PARAMETER :: KARMAN = 0.4\r
1246 \r
1247 !... Integer\r
1248 INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L\r
1249 INTEGER :: KCBLMX\r
1250 INTEGER, DIMENSION( its:ite ) :: KCBL\r
1251 \r
1252 !... Real\r
1253 REAL :: MBMAX, HOVL, MEDDY, MBAR\r
1254 REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1\r
1255 REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, DTLIM\r
1256 REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN\r
1257 REAL, DIMENSION( 1:NSP, its:ite ) :: FS\r
1258 REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS\r
1259 REAL DELC\r
1260 REAL, DIMENSION( 1:NSP,its:ite,kts:kte ) :: VCI\r
1261 \r
1262 REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI !, Y\r
1263 REAL, DIMENSION( 1:NSP,kts:kte ) :: DI, UI \r
1264 !\r
1265 !--Start Exicutable ----\r
1266 \r
1267 ILX = ite\r
1268 KL = kte\r
1269 KLM = kte - 1\r
1270 \r
1271 KCBLMX = 0\r
1272 MBMAX = 0.0\r
1273 \r
1274 !---COMPUTE ACM MIXING RATE\r
1275 DO I = its, ILX\r
1276 DTLIM(I) = DTPBL\r
1277 PSTARI(I) = 1.0 / PSTAR(I)\r
1278 KCBL(I) = 1\r
1279 FSACM(I) = 0.0\r
1280 \r
1281 IF (NOCONV(I) .EQ. 1) THEN\r
1282 KCBL(I) = KLPBL(I)\r
1283 \r
1284 !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE\r
1285 !--New couple ACM & EDDY-------------------------------------------------------------\r
1286 HOVL = -PBL(I) / MOL(I)\r
1287 FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))\r
1288 MEDDY = EDDYZM(I,1) * DZFI(i,1) / (PBL(I) - ZF(i,1))\r
1289 MBAR = MEDDY * FSACM(I)\r
1290 DO K = kts,KCBL(I)-1\r
1291 EDDYZM(I,K) = EDDYZM(I,K) * (1.0 - FSACM(I))\r
1292 ENDDO\r
1293 \r
1294 MBMAX = AMAX1(MBMAX,MBAR)\r
1295 DO K = kts+1,KCBL(I)\r
1296 MBARKS(K,I) = MBAR\r
1297 MDWN(K,I) = MBAR * (PBL(I) - ZF(i,K-1)) * DZHI(i,K)\r
1298 ENDDO\r
1299 MBARKS(1,I) = MBAR\r
1300 MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)\r
1301 MDWN(KCBL(I)+1,I) = 0.0\r
1302 ENDIF\r
1303 ENDDO ! end of I loop\r
1304 \r
1305 DO K = kts,KLM\r
1306 DO I = its,ILX\r
1307 EKZ = EDDYZM(I,K) * DZFI(i,K) * DZHI(i,K)\r
1308 DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))\r
1309 ENDDO\r
1310 ENDDO\r
1311 \r
1312 DO I = its,ILX \r
1313 IF (NOCONV(I) .EQ. 1) THEN\r
1314 KCBLMX = AMAX0(KLPBL(I),KCBLMX)\r
1315 RZ = (ZF(i,KCBL(I)) - ZF(i,1)) * DZHI(i,1)\r
1316 DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))\r
1317 ENDIF\r
1318 ENDDO\r
1319 \r
1320 DO K = kts,KL\r
1321 DO I = its,ILX\r
1322 VCI(1,I,K) = US(I,K)\r
1323 VCI(2,I,K) = VS(I,K)\r
1324 ENDDO\r
1325 ENDDO\r
1326 \r
1327 NSPX=2\r
1328 \r
1329 DO I = its,ILX\r
1330 FM = -USTM(I) * USTM(I)\r
1331 WSPD = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9\r
1332 FS(1,I) = FM * US(I,1) / WSPD\r
1333 FS(2,I) = FM * VS(I,1) / WSPD\r
1334 ENDDO\r
1335 !\r
1336 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
1337 DO I = its,ILX \r
1338 \r
1339 NLP = INT(DTPBL / DTLIM(I) + 1.0)\r
1340 DTS = (DTPBL / NLP)\r
1341 DTRAT = DTS / DTPBL\r
1342 DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP \r
1343 \r
1344 !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES\r
1345 \r
1346 DO K = kts,KL\r
1347 AI(K) = 0.0\r
1348 BI(K) = 0.0\r
1349 CI(K) = 0.0\r
1350 EI(K) = 0.0\r
1351 ENDDO\r
1352 \r
1353 DO K = 2, KCBL(I)\r
1354 EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DZH(i,K) * DZHI(i,K-1)\r
1355 BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS\r
1356 AI(K) = -CRANKP * MBARKS(K,I) * DTS\r
1357 ENDDO\r
1358 \r
1359 EI(1) = EI(1) -EDDYZM(I,1) * CRANKP * DZHI(i,1) * DZFI(i,1) * DTS\r
1360 AI(2) = AI(2) -EDDYZM(I,1) * CRANKP * DZHI(i,2) * DZFI(i,1) * DTS\r
1361 \r
1362 DO K = KCBL(I)+1, KL\r
1363 BI(K) = 1.0\r
1364 ENDDO\r
1365 \r
1366 DO K = 2,KL\r
1367 XPLUS(K) = EDDYZM(I,K) * DZHI(i,K) * DZFI(i,K) * DTS\r
1368 XMINUS(K) = EDDYZM(I,K-1) * DZHI(i,K) * DZFI(i,K-1) * DTS\r
1369 CI(K) = - XMINUS(K) * CRANKP\r
1370 EI(K) = EI(K) - XPLUS(K) * CRANKP\r
1371 BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP\r
1372 ENDDO\r
1373 \r
1374 IF (NOCONV(I) .EQ. 1) THEN\r
1375 BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBL(I) - ZF(i,1)) * DTS &\r
1376 * DZHI(i,1) + EDDYZM(I,1) * DZHI(i,1) * DZFI(i,1) * CRANKP * DTS\r
1377 ELSE\r
1378 BI(1) = 1.0 + EDDYZM(I,1) * DZHI(i,1) * DZFI(i,1) * CRANKP * DTS\r
1379 ENDIF\r
1380 \r
1381 \r
1382 DO K = 1,KL\r
1383 DO L = 1,NSPX \r
1384 DI(L,K) = 0.0\r
1385 ENDDO\r
1386 ENDDO\r
1387 !\r
1388 !** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION\r
1389 DO K = 2,KCBL(I)\r
1390 DO L = 1,NSPX \r
1391 DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * &\r
1392 VCI(L,I,K) + DZH(i,K+1) * DZHI(i,K) * &\r
1393 MDWN(K+1,I) * VCI(L,I,K+1))\r
1394 DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC\r
1395 ENDDO\r
1396 ENDDO\r
1397 \r
1398 DO K = KCBL(I)+1, KL\r
1399 DO L = 1,NSPX \r
1400 DI(L,K) = VCI(L,I,K)\r
1401 ENDDO\r
1402 ENDDO\r
1403 \r
1404 DO K = 2,KL\r
1405 IF (K .EQ. KL) THEN\r
1406 DO L = 1,NSPX \r
1407 DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * &\r
1408 (VCI(L,I,K) - VCI(L,I,K-1))\r
1409 ENDDO\r
1410 ELSE\r
1411 DO L = 1,NSPX \r
1412 DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * &\r
1413 (VCI(L,I,K+1) - VCI(L,I,K)) - &\r
1414 (1.0 - CRANKP) * XMINUS(K) * &\r
1415 (VCI(L,I,K) - VCI(L,I,K-1))\r
1416 ENDDO\r
1417 ENDIF\r
1418 ENDDO\r
1419 \r
1420 IF (NOCONV(I) .EQ. 1) THEN\r
1421 DO L = 1,NSPX \r
1422 DI(L,1) = VCI(L,I,1) + (FS(L,I) - (1.0 - CRANKP) &\r
1423 * (MBARKS(1,I) * &\r
1424 (PBL(I) - ZF(i,1)) * VCI(L,I,1) - &\r
1425 MDWN(2,I) * VCI(L,I,2) * DZH(i,2))) * DZHI(i,1) * DTS\r
1426 ENDDO\r
1427 ELSE\r
1428 DO L = 1,NSPX \r
1429 DI(L,1) = VCI(L,I,1) + FS(L,I) * DZHI(i,1) * DTS\r
1430 ENDDO\r
1431 ENDIF\r
1432 DO L = 1,NSPX \r
1433 DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZM(I,1) * DZHI(i,1) &\r
1434 * DZFI(i,1) * DTS * (VCI(L,I,2) - VCI(L,I,1))\r
1435 ENDDO\r
1436 IF ( NOCONV(I) .EQ. 1 ) THEN\r
1437 CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)\r
1438 ELSE\r
1439 CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)\r
1440 END IF\r
1441 !\r
1442 !-- COMPUTE NEW THETAV AND Q\r
1443 DO K = 1,KL\r
1444 DO L = 1,NSPX \r
1445 VCI(L,I,K) = UI(L,K)\r
1446 ENDDO\r
1447 ENDDO\r
1448 \r
1449 ENDDO ! END I LOOP\r
1450 ENDDO ! END SUB TIME LOOP\r
1451 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
1452 \r
1453 !\r
1454 DO K = kts,KL\r
1455 DO I = its,ILX\r
1456 UX(I,K) = VCI(1,I,K)\r
1457 VX(I,K) = VCI(2,I,K)\r
1458 ENDDO\r
1459 ENDDO\r
1460 \r
1461 END SUBROUTINE ACMM\r
1462 !-----------------------------------------------------------------------\r
1463 \r
1464 !-----------------------------------------------------------------------\r
1465 !-----------------------------------------------------------------------\r
1466 SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)\r
1467 \r
1468 !-----------------------------------------------------------------------\r
1469 !-----------------------------------------------------------------------\r
1470 IMPLICIT NONE\r
1471 !\r
1472 !-- Bordered band diagonal matrix solver for ACM2\r
1473 \r
1474 !-- ACM2 Matrix is in this form:\r
1475 ! B1 E1\r
1476 ! A2 B2 E2\r
1477 ! A3 C3 B3 E3\r
1478 ! A4 C4 B4 E4\r
1479 ! A5 C5 B5 E5\r
1480 ! A6 C6 B6\r
1481 \r
1482 !--Upper Matrix is\r
1483 ! U11 U12\r
1484 ! U22 U23\r
1485 ! U33 U34\r
1486 ! U44 U45\r
1487 ! U55 U56\r
1488 ! U66\r
1489 \r
1490 !--Lower Matrix is:\r
1491 ! 1\r
1492 ! L21 1\r
1493 ! L31 L32 1\r
1494 ! L41 L42 L43 1\r
1495 ! L51 L52 L53 L54 1\r
1496 ! L61 L62 L63 L64 L65 1\r
1497 !---------------------------------------------------------\r
1498 !...Arguments\r
1499 INTEGER, INTENT(IN) :: KL\r
1500 INTEGER, INTENT(IN) :: NSP\r
1501 REAL A(KL),B(KL),E(KL)\r
1502 REAL C(KL),D(NSP,KL),X(NSP,KL)\r
1503 \r
1504 !...Locals\r
1505 REAL Y(NSP,KL),AIJ,SUM\r
1506 REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)\r
1507 INTEGER I,J,V\r
1508 \r
1509 !-- Define Upper and Lower matrices\r
1510 L(1,1) = 1.\r
1511 UII(1) = B(1)\r
1512 RUII(1) = 1./UII(1)\r
1513 DO I = 2, KL\r
1514 L(I,I) = 1.\r
1515 L(I,1) = A(I)/B(1)\r
1516 UIIP1(I-1)=E(I-1)\r
1517 IF(I.GE.3) THEN\r
1518 DO J = 2,I-1\r
1519 IF(I.EQ.J+1) THEN\r
1520 AIJ = C(I)\r
1521 ELSE\r
1522 AIJ = 0.\r
1523 ENDIF\r
1524 L(I,J) = (AIJ-L(I,J-1)*E(J-1))/ &\r
1525 (B(J)-L(J,J-1)*E(J-1))\r
1526 ENDDO\r
1527 ENDIF\r
1528 ENDDO\r
1529 \r
1530 DO I = 2,KL\r
1531 UII(I) = B(I)-L(I,I-1)*E(I-1)\r
1532 RUII(I) = 1./UII(I)\r
1533 ENDDO\r
1534 \r
1535 !-- Forward sub for Ly=d\r
1536 DO V= 1, NSP\r
1537 Y(V,1) = D(V,1)\r
1538 DO I=2,KL\r
1539 SUM = D(V,I)\r
1540 DO J=1,I-1\r
1541 SUM = SUM - L(I,J)*Y(V,J)\r
1542 ENDDO\r
1543 Y(V,I) = SUM\r
1544 ENDDO\r
1545 ENDDO\r
1546 \r
1547 !-- Back sub for Ux=y\r
1548 \r
1549 DO V= 1, NSP\r
1550 X(V,KL) = Y(V,KL)*RUII(KL)\r
1551 ENDDO\r
1552 DO I = KL-1,1,-1\r
1553 DO V= 1, NSP\r
1554 X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I)\r
1555 ENDDO\r
1556 ENDDO\r
1557 \r
1558 END SUBROUTINE MATRIX\r
1559 \r
1560 \r
1561 !-----------------------------------------------------------------------\r
1562 !-----------------------------------------------------------------------\r
1563 SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)\r
1564 !-----------------------------------------------------------------------\r
1565 !-----------------------------------------------------------------------\r
1566 \r
1567 ! FUNCTION:\r
1568 ! Solves tridiagonal system by Thomas algorithm. \r
1569 ! The associated tri-diagonal system is stored in 3 arrays\r
1570 ! D : diagonal\r
1571 ! L : sub-diagonal\r
1572 ! U : super-diagonal\r
1573 ! B : right hand side function\r
1574 ! X : return solution from tridiagonal solver\r
1575 \r
1576 ! [ D(1) U(1) 0 0 0 ... 0 ]\r
1577 ! [ L(2) D(2) U(2) 0 0 ... . ]\r
1578 ! [ 0 L(3) D(3) U(3) 0 ... . ]\r
1579 ! [ . . . . . ] X(i) = B(i)\r
1580 ! [ . . . . 0 ]\r
1581 ! [ . . . . ]\r
1582 ! [ 0 L(n) D(n) ]\r
1583 \r
1584 !-----------------------------------------------------------------------\r
1585 \r
1586 IMPLICIT NONE\r
1587 \r
1588 ! Arguments:\r
1589 \r
1590 INTEGER, INTENT(IN) :: KL\r
1591 INTEGER, INTENT(IN) :: NSP\r
1592 \r
1593 REAL L( KL ) ! subdiagonal\r
1594 REAL D(KL) ! diagonal\r
1595 REAL U( KL ) ! superdiagonal\r
1596 REAL B(NSP,KL ) ! R.H. side\r
1597 REAL X( NSP,KL ) ! solution\r
1598 \r
1599 ! Local Variables:\r
1600 \r
1601 REAL GAM( KL )\r
1602 REAL BET\r
1603 INTEGER V, K\r
1604 \r
1605 ! Decomposition and forward substitution:\r
1606 BET = 1.0 / D( 1 )\r
1607 DO V = 1, NSP\r
1608 X( V,1 ) = BET * B(V,1 )\r
1609 ENDDO\r
1610 \r
1611 DO K = 2, KL\r
1612 GAM(K ) = BET * U( K-1 )\r
1613 BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )\r
1614 DO V = 1, NSP\r
1615 X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )\r
1616 ENDDO\r
1617 ENDDO\r
1618 \r
1619 ! Back-substitution:\r
1620 \r
1621 DO K = KL - 1, 1, -1\r
1622 DO V = 1, NSP\r
1623 X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )\r
1624 ENDDO\r
1625 ENDDO\r
1626 \r
1627 END SUBROUTINE TRI\r
1628 !-----------------------------------------------------------------------\r
1629 !-----------------------------------------------------------------------\r
1630 \r
1631 END MODULE module_bl_acm\r
1632 \r