Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_acm.F @ 3431

Last change on this file since 3431 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 43.4 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE module_bl_acm
4
5	! USE module_model_constants
6
7	REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number
8	REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER
9
10	CONTAINS
11
12	!-------------------------------------------------------------------
13	SUBROUTINE ACMPBL(XTIME, DTPBL, ZNW, SIGMAH, &
14	U3D, V3D, PP3D, DZ8W, TH3D, T3D, &
15	QV3D, QC3D, QI3D, RR3D, & !! Moisture
16	UST, HFX, QFX, TSK, &
17	PSFC, EP1, G, &
18	ROVCP, RD, CPD, &
19	PBLH, KPBL2D, REGIME, &
20	GZ1OZ0, WSPD, PSIM, MUT, &
21	RUBLTEN, RVBLTEN, RTHBLTEN, &
22	RQVBLTEN, RQCBLTEN, RQIBLTEN, &
23	ids,ide, jds,jde, kds,kde, &
24	ims,ime, jms,jme, kms,kme, &
25	its,ite, jts,jte, kts,kte)
26	!**********************************************************************
27
28	! THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO
29	! THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2 (ACM2), WHICH IS A COMBINED
30	! LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
31	!
32	! REFERENCES:
33	! Pleim (2007) A combined local and non-local closure model for the atmospheric
34	! boundary layer. Part1: Model description and testing.
35	! JAMC, 46, 1383-1395
36	! Pleim (2007) A combined local and non-local closure model for the atmospheric
37	! boundary layer. Part2: Application and evaluation in a mesoscale
38	! meteorology model. JAMC, 46, 1396-1409
39	!
40	! REVISION HISTORY:
41	! AX 3/2005 - developed WRF version based on the MM5 PX LSM
42	! RG and JP 7/2006 - Finished WRF adaptation
43	!
44	!**********************************************************************
45	! ARGUMENT LIST:
46	!
47	!... Inputs:
48	!-- XTIME Time since simulation start (min)
49	!-- DTPBL PBL time step
50	!-- ZNW Sigma at full layer
51	!-- SIGMAH Sigma at half layer
52	!-- U3D 3D u-velocity interpolated to theta points (m/s)
53	!-- V3D 3D v-velocity interpolated to theta points (m/s)
54	!-- PP3D Pressure at half levels (Pa)
55	!-- DZ8W dz between full levels (m)
56	!-- TH3D Potential Temperature (K)
57	!-- T3D Temperature (K)
58	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
59	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
60	!-- QI3D 3D ice mixing ratio (Kg/Kg)
61	!-- RR3D 3D dry air density (kg/m^3)
62	!-- UST Friction Velocity (m/s)
63	!-- HFX Upward heat flux at the surface (w/m^2)
64	!-- QFX Upward moisture flux at the surface (Kg/m^2/s)
65	!-- TSK Surface temperature (K)
66	!-- PSFC Pressure at the surface (Pa)
67	!-- EP1 Constant for virtual temperature (r_v/r_d-1) (dimensionless)
68	!-- G Gravity (m/s^2)
69	!-- ROVCP r/cp
70	!-- RD gas constant for dry air (j/kg/k)
71	!-- CPD heat capacity at constant pressure for dry air (j/kg/k)
72	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
73	!-- WSPD wind speed at lowest model level (m/s)
74	!-- PSIM similarity stability function for momentum
75	!-- MUT Total Mu : Psfc - Ptop
76
77	!-- ids start index for i in domain
78	!-- ide end index for i in domain
79	!-- jds start index for j in domain
80	!-- jde end index for j in domain
81	!-- kds start index for k in domain
82	!-- kde end index for k in domain
83	!-- ims start index for i in memory
84	!-- ime end index for i in memory
85	!-- jms start index for j in memory
86	!-- jme end index for j in memory
87	!-- kms start index for k in memory
88	!-- kme end index for k in memory
89	!-- jts start index for j in tile
90	!-- jte end index for j in tile
91	!-- kts start index for k in tile
92	!-- kte end index for k in tile
93	!
94	!... Outputs:
95	!-- PBLH PBL height (m)
96	!-- KPBL2D K index for PBL layer
97	!-- REGIME Flag indicating PBL regime (stable, unstable, etc.)
98	!-- RUBLTEN U tendency due to PBL parameterization (m/s^2)
99	!-- RVBLTEN V tendency due to PBL parameterization (m/s^2)
100	!-- RTHBLTEN Theta tendency due to PBL parameterization (K/s)
101	!-- RQVBLTEN Qv tendency due to PBL parameterization (kg/kg/s)
102	!-- RQCBLTEN Qc tendency due to PBL parameterization (kg/kg/s)
103	!-- RQIBLTEN Qi tendency due to PBL parameterization (kg/kg/s)
104	!
105	!-------------------------------------------------------------------
106	IMPLICIT NONE
107
108	!.......Arguments
109	! DECLARATIONS - INTEGER
110	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
111	ims,ime, jms,jme, kms,kme, &
112	its,ite, jts,jte, kts,kte, XTIME
113
114	! DECLARATIONS - REAL
115	REAL, INTENT(IN) :: DTPBL, EP1, &
116	G, ROVCP, RD, CPD
117	REAL, DIMENSION( kms:kme ), INTENT(IN) :: ZNW, SIGMAH
118	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
119	INTENT(IN) :: U3D, V3D, &
120	PP3D, DZ8W, T3D, &
121	QV3D, QC3D, QI3D, &
122	RR3D, TH3D
123	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
124	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
125	HFX, QFX, TSK, &
126	PSFC
127	INTEGER, DIMENSION( ims:ime, jms:jme ), &
128	INTENT(OUT ) :: KPBL2D
129	REAL, DIMENSION( ims:ime, jms:jme ), &
130	INTENT(INOUT) :: PBLH
131	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: REGIME
132	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
133	psim, &
134	gz1oz0, &
135	wspd , mut
136	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
137	INTENT(INOUT) :: RUBLTEN, RVBLTEN, &
138	RTHBLTEN, RQVBLTEN, &
139	RQCBLTEN, RQIBLTEN
140
141	!... Local Variables
142
143	!... Integer
144	INTEGER :: J, K
145	!... Real
146	REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
147	REAL, DIMENSION( 0:kte ) :: SIGMAF
148	REAL RDT
149	REAL, PARAMETER :: KARMAN = 0.4
150	!...
151
152	RDT = 1.0 / DTPBL
153
154	DO K = 1, kte
155	SIGMAF(K-1) = ZNW(K)
156	ENDDO
157	SIGMAF(kte) = 0.0
158
159	DO K = 1, kte
160	DSIGH(K) = SIGMAF(K) - SIGMAF(K-1)
161	DSIGHI(K) = 1.0 / DSIGH(K)
162	ENDDO
163
164	DO K = kts,kte-1
165	DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
166	ENDDO
167
168	DSIGFI(kte) = DSIGFI(kte-1)
169	do j = jts,jte
170	CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH &
171	,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH &
172	,us=u3d(ims,kms,j),vs=v3d(ims,kms,j) &
173	,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j) &
174	,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j) &
175	,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j) &
176	,densx=RR3D(ims,kms,j) &
177	,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
178	,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j) &
179	,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
180	,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt &
181	,psfcpa=psfc(ims,j),ust=ust(ims,j) &
182	,pbl=pblh(ims,j) &
183	,regime=regime(ims,j),psim=psim(ims,j) &
184	,hfx=hfx(ims,j),qfx=qfx(ims,j) &
185	,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &
186	,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j) &
187	,mut=mut(ims,j) &
188	,ep1=ep1,karman=karman &
189	,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
190	,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
191	,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
192	enddo
193
194	END SUBROUTINE ACMPBL
195
196	!--------------------------------------------------------------------
197	SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah &
198	,dsigfi,dsighi,dsigh &
199	,us,vs,theta,tt,qvs,qcs,qis &
200	,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp &
201	,cpd,g,rovcp,rd,rdt,psfcpa,ust &
202	,pbl,regime,psim &
203	,hfx,qfx,tg,gz1oz0,wspd ,klpbl &
204	,mut &
205	,ep1,karman &
206	,ids,ide, jds,jde, kds,kde &
207	,ims,ime, jms,jme, kms,kme &
208	,its,ite, jts,jte, kts,kte )
209
210	IMPLICIT NONE
211
212	!.......Arguments
213
214	!... Real
215	REAL, DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
216	REAL, DIMENSION( kms:kme ), INTENT(IN) :: SIGMAH
217	REAL, DIMENSION( kts:kte ), INTENT(IN) :: DSIGH, DSIGHI, DSIGFI
218	REAL , INTENT(IN) :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
219	REAL , DIMENSION( ims:ime ), INTENT(INOUT) :: PBL, UST
220
221	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, TT, &
222	QVS, QCS, QIS, DENSX
223	REAL, DIMENSION( ims:ime, kms:kme ), intent(in) :: DZF
224	REAL, DIMENSION( ims:ime, kms:kme ), intent(inout) :: utnp, &
225	vtnp, &
226	ttnp, &
227	qvtnp, &
228	qctnp, &
229	qitnp
230	real, dimension( ims:ime ), intent(in ) :: psfcpa
231	real, dimension( ims:ime ), intent(in ) :: tg
232	real, dimension( ims:ime ), intent(inout) :: regime
233	real, dimension( ims:ime ), intent(in) :: wspd, psim, gz1oz0
234	real, dimension( ims:ime ), intent(in) :: hfx, qfx
235	real, dimension( ims:ime ), intent(in) :: mut
236
237	!... Integer
238	INTEGER, DIMENSION( ims:ime ), INTENT(OUT):: KLPBL
239	INTEGER, INTENT(IN) :: XTIME
240	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
241	ims,ime, jms,jme, kms,kme, &
242	its,ite, jts,jte, kts,kte, j
243	!--------------------------------------------------------------------
244	!--Local
245	INTEGER I, K
246	INTEGER :: KPBLHT
247	INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV
248
249	!... Real
250	REAL :: TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
251	REAL :: psix, THV1
252	REAL, DIMENSION( its:ite ) :: FINT, PSTAR, CPAIR
253	REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB, &
254	EDDYZ, UX, VX, THETAX, &
255	QVX, QCX, QIX, ZA
256	REAL, DIMENSION( its:ite, 0:kte ) :: ZF
257	REAL, DIMENSION( its:ite) :: WST, TST, QST, USTM, TSTV
258	REAL, DIMENSION( its:ite ) :: PBLSIG, MOL
259	REAL :: FINTT, ZMIX, UMIX, VMIX
260	REAL :: TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
261	REAL :: A,TST12,RL,ZFUNC
262	! REAL, PARAMETER :: KARMAN = 0.4
263
264	!... Integer
265	INTEGER :: KL, jtf, ktf, itf, KMIX
266	!...
267	!-----initialize vertical tendencies and
268	!
269	do i = its,ite
270	do k = kts,kte
271	utnp(i,k) = 0.
272	vtnp(i,k) = 0.
273	ttnp(i,k) = 0.
274	enddo
275	enddo
276	!
277	do k = kts,kte
278	do i = its,ite
279	qvtnp(i,k) = 0.
280	enddo
281	enddo
282	!
283	do k = kts,kte
284	do i = its,ite
285	qctnp(i,k) = 0.
286	qitnp(i,k) = 0.
287	enddo
288	enddo
289	!
290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291	!!! Compute Micromet Scaling variables, not availiable in WRF for ACM
292	DO I = its,ite
293	cpair(i) = CPD * (1.0 + 0.84 * QVS(I,1)) ! J/(K KG)
294	TMPFX = HFX(I) / (cpair(i) * DENSX(I,1))
295	TMPVTCON = 1.0 + EP1 * QVS(I,1) ! COnversion factor for virtual temperature
296	WS1 = SQRT(US(I,1)2 + VS(I,1)2) ! Level 1 wind speed
297	TST(I) = -TMPFX / UST(I)
298	QST(I) = -QFX(I) / (UST(I)*DENSX(I,1))
299	USTM(I) = UST(I) * WS1 / wspd(i)
300	THV1 = TMPVTCON * THETA(I,1)
301	TSTV(I) = TST(I)TMPVTCON + THV1EP1*QST(I)
302	MOL(I) = THV1 * UST(i)*2/(KARMANG*TSTV(I))
303	WST(I) = UST(I) * (PBL(I)/(KARMANABS(MOL(I)))) * 0.333333
304	PSTAR(I) = MUT(I)/1000. ! P* in cb
305	ENDDO
306	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307
308	!... Compute PBL height
309
310	!... compute the height of full- and half-sigma level above ground level
311	DO I = its,ite
312	ZF(I,0) = 0.0
313	KLPBL(I) = 1
314	ENDDO
315
316	DO K = kts, kte
317	DO I = its,ite
318	ZF(I,K) = DZF(I,K) + ZF(I,K-1)
319	ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
320	ENDDO
321	ENDDO
322
323	DO K = kts, kte
324	DO I = its,ite
325	TVCON = 1.0 + EP1 * QVS(I,K)
326	THETAV(I,K) = THETA(I,K) * TVCON
327	ENDDO
328	ENDDO
329
330
331	!... COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990
332	DO 100 I = its,ite
333	if(MOL(I).LT.0.0) then
334	WSS = (UST(I) ** 3 + 0.6 * WST(I) 3) 0.33333
335	TMPFX = -UST(I) * TSTV(I)
336	TCONV = 8.5 * TMPFX / WSS
337	TH1 = THETAV(I,1) + TCONV
338	else
339	TH1 = THETAV(I,1)
340	endif
341	99 KMIX = 1
342	DO K = 1,kte
343	DTMP = THETAV(I,K) - TH1
344	IF (DTMP.LT.0.0) KMIX = K
345	ENDDO
346	IF(KMIX.GT.1) THEN
347	FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1) &
348	- THETAV(I,KMIX))
349	ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
350	UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
351	VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
352	ELSE
353	ZMIX = ZA(I,1)
354	UMIX = US(I,1)
355	VMIX = VS(I,1)
356	ENDIF
357	DO K = KMIX,kte
358	DTMP = THETAV(I,K) - TH1
359	TOG = 0.5 * (THETAV(I,K) + TH1) / G
360	WSSQ = (US(I,K)-UMIX)**2 &
361	+ (VS(I,K)-VMIX)**2
362	IF (KMIX == 1) WSSQ = WSSQ + 100.UST(I)UST(I)
363	WSSQ = MAX( WSSQ, 0.1 )
364	RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
365	IF (RIB(I,K) .GE. RIC) GO TO 201
366	ENDDO
367
368	print *,' RIB never exceeds RIC, RIB(i,kte) = ',rib(i,5), &
369	' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i), &
370	' TCONV = ',TCONV,' WST = ',WST(I), &
371	' KMIX = ',kmix,' UST = ',UST(I), &
372	' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1), &
373	' I,J=',I,J
374	STOP
375	201 CONTINUE
376
377	KPBLH(I) = K
378
379	100 CONTINUE
380
381	DO I = its,ite
382	IF (KPBLH(I) .NE. 1) THEN
383	!---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
384	FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) - &
385	RIB(I,KPBLH(I)-1))
386	IF (FINT(I) .GT. 0.5) THEN
387	KPBLHT = KPBLH(I)
388	FINT(I) = FINT(I) - 0.5
389	ELSE
390	KPBLHT = KPBLH(I) - 1
391	FINT(I) = FINT(I) + 0.5
392	ENDIF
393	PBL(I) = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) + &
394	ZF(I,KPBLHT-1)
395	KLPBL(I) = KPBLHT
396	PBLSIG(I) = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1) ! sigma at PBL height
397	ELSE
398	KLPBL(I) = 1
399	PBL(I) = ZF(I,1)
400	PBLSIG(I) = SIGMAF(1)
401	ENDIF
402
403	ENDDO
404
405	DO I = its,ite
406	NOCONV(I) = 0
407
408	! Check for CBL and identify conv. vs. non-conv cells
409	IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3 &
410	.AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
411	NOCONV(I) = 1
412	REGIME(I) = 4.0 ! FREE CONVECTIVE - ACM
413	ENDIF
414	ENDDO
415
416	!... Calculate Kz
417	CALL EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
418	US, VS, TT, THETAV, DENSX, PSTAR, &
419	QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
420	EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)
421
422	CALL ACM(DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, J, &
423	KLPBL, PBL, PBLSIG, MOL, UST, &
424	TST, QST, USTM, EDDYZ, DENSX, &
425	US, VS, THETA, QVS, QCS, QIS, &
426	UX, VX, THETAX, QVX, QCX, QIX, &
427	ids,ide, jds,jde, kds,kde, &
428	ims,ime, jms,jme, kms,kme, &
429	its,ite, jts,jte, kts,kte)
430
431	!... Calculate tendency due to PBL parameterization
432
433	DO K = kts, kte
434	DO I = its, ite
435	UTNP(I,K) = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
436	VTNP(I,K) = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
437	TTNP(I,K) = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
438	QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
439	QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
440	QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
441	ENDDO
442	ENDDO
443	! IF(J.EQ.36)PRINT *,' PBL,THETA,THETAX,QVS,QVX=',PBL(20),THETA(20,1),THETAX(20,1),QVS(20,1),QVX(20,1)
444	! IF(J.EQ.36)PRINT *,' UST,ustm,TG,TST=',UST(20),ustm(20),Tg(20),TST(20)
445	! IF(J.EQ.36)PRINT *,' HFX,MOL,WST,WSPD=',HFX(20),MOL(20),WST(20),wspd(20)
446	! IF(J.EQ.36)then
447	! i=20
448	! do k=1,kte
449	! PRINT *,' qvten,uten,vten=',QVTNP(I,K),UTNP(I,K),VTNP(I,K)
450	! PRINT *,' k,thten,th,thx,edy=',k,TTNP(I,K),THETA(20,k),THETAX(20,k),EDDYZ(20,K)
451	! enddo
452	! ENDIF
453	END SUBROUTINE ACM2D
454
455	!-------------------------------------------------------------------
456	SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
457	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
458	restart, allowed_to_read , &
459	ids, ide, jds, jde, kds, kde, &
460	ims, ime, jms, jme, kms, kme, &
461	its, ite, jts, jte, kts, kte )
462	!**********************************************************************
463	!
464	! This subroutine is for preparing ACM PBL variables.
465	! Called from module_physics_init.F
466	!
467	! REVISION HISTORY:
468	! AX 3/2005 - Originally developed
469	!**********************************************************************
470	! ARGUMENT LIST:
471	!
472	!-------------------------------------------------------------------
473	IMPLICIT NONE
474	!
475	LOGICAL , INTENT(IN) :: restart , allowed_to_read
476
477	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
478	ims, ime, jms, jme, kms, kme, &
479	its, ite, jts, jte, kts, kte
480
481	INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
482
483	! REAL , DIMENSION( kms:kme ), INTENT(IN) :: SHALF
484	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
485	RUBLTEN, &
486	RVBLTEN, &
487	RTHBLTEN, &
488	RQVBLTEN, &
489	RQCBLTEN, &
490	RQIBLTEN
491
492	!... Local Variables
493	INTEGER :: i, j, k, itf, jtf, ktf
494
495	!
496	jtf=min0(jte,jde-1)
497	ktf=min0(kte,kde-1)
498	itf=min0(ite,ide-1)
499
500	IF(.not.restart)THEN
501	DO j=jts,jtf
502	DO k=kts,ktf
503	DO i=its,itf
504	RUBLTEN(i,k,j)=0.
505	RVBLTEN(i,k,j)=0.
506	RTHBLTEN(i,k,j)=0.
507	RQVBLTEN(i,k,j)=0.
508	RQCBLTEN(i,k,j)=0.
509	ENDDO
510	ENDDO
511	ENDDO
512	ENDIF
513
514	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
515	DO j=jts,jtf
516	DO k=kts,ktf
517	DO i=its,itf
518	RQIBLTEN(i,k,j)=0.
519	ENDDO
520	ENDDO
521	ENDDO
522	ENDIF
523
524
525	END SUBROUTINE acminit
526
527	!-------------------------------------------------------------------
528	SUBROUTINE EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
529	US, VS, TT, THETAV, DENSX, PSTAR, &
530	QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
531	EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )
532
533
534	!**********************************************************************
535	! Two methods for computing Kz:
536	! 1. Boundary scaling similar to Holtslag and Boville (1993)
537	! 2. Local Kz computed as function of local Richardson # and vertical
538	! wind shear, similar to LIU & CARROLL (1996)
539	!
540	!**********************************************************************
541	!
542	!-- DTPBL time step of the minor loop for the land-surface/pbl model
543	!-- ZF height of full sigma level
544	!-- ZA height of half sigma level
545	!-- MOL Monin-Obukhov length in 1D form
546	!-- PBL PBL height in 1D form
547	!-- UST friction velocity U* in 1D form (m/s)
548	!-- US U wind
549	!-- VS V wind
550	!-- TT temperature
551	!-- THETAV potential virtual temperature
552	!-- DENSX dry air density (kg/m^3)
553	!-- PSTAR P*=Psfc-Ptop
554	!-- QVS water vapor mixing ratio (Kg/Kg)
555	!-- QCS cloud mixing ratio (Kg/Kg)
556	!-- QIS ice mixing ratio (Kg/Kg)
557	!-- DSIGFI inverse of sigma layer delta
558	!-- G gravity
559	!-- RD gas constant for dry air (j/kg/k)
560	!-- CPAIR specific heat of moist air (M^2 S^-2 K^-1)
561	!-- EDDYZ eddy diffusivity KZ
562	!-----------------------------------------------------------------------
563
564	IMPLICIT NONE
565
566	!.......Arguments
567
568	!... Integer
569	INTEGER, INTENT(IN) :: its,ite, kts,kte,ims,ime,kms,kme
570	!... Real
571	REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
572	REAL , INTENT(IN) :: DTPBL, G, RD
573	REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGFI
574	REAL , DIMENSION( its:ite ), INTENT(IN) :: MOL, PSTAR, CPAIR
575
576	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, TT, &
577	QVS, QCS, QIS, DENSX
578	REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
579	REAL, DIMENSION( its:ite, 0:kte ) , INTENT(IN) :: ZF
580
581	REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ
582
583	!.......Local variables
584
585	!... Integer
586	INTEGER :: ILX, KL, KLM, K, I
587
588	!... Real
589	REAL :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
590	REAL :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF,kzo
591
592	!... Parameters
593	REAL, PARAMETER :: RV = 461.5
594	REAL, PARAMETER :: RC = 0.25
595	REAL, PARAMETER :: RLAM = 80.0
596	REAL, PARAMETER :: GAMH = 16.0 !15.0 ! Holtslag and Boville (1993)
597	REAL, PARAMETER :: BETAH = 5.0 ! Holtslag and Boville (1993)
598	REAL, PARAMETER :: KARMAN = 0.4
599	REAL, PARAMETER :: EDYZ0 = 0.01 ! New Min Kz
600	! REAL, PARAMETER :: EDYZ0 = 0.1
601	!-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
602	INTEGER, PARAMETER :: imvdif = 1
603	!
604	ILX = ite
605	KL = kte
606	KLM = kte - 1
607
608	DO K = kts,KLM
609	DO I = its,ILX
610	EDYZ = 0.0
611	ZOVL = 0.0
612	DZF = ZA(I,K+1) - ZA(I,K)
613	! kzo = min(0.001 * dzf,EDYZ0)
614	kzo = 0.001 * dzf
615	! kzo = EDYZ0
616	!-------------------------------------------------
617	IF (ZF(I,K) .LT. PBL(I)) THEN
618	ZOVL = ZF(I,K) / MOL(I)
619	IF (ZOVL .LT. 0.0) THEN
620	IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
621	PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
622	WT = UST(I) / PHIH
623	ELSE
624	ZSOL = 0.1 * PBL(I) / MOL(I)
625	PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
626	WT = UST(I) / PHIH
627	ENDIF
628	ELSE IF (ZOVL .LT. 1.0) THEN
629	PHIH = 1.0 + BETAH * ZOVL
630	WT = UST(I) / PHIH
631	ELSE
632	PHIH = BETAH + ZOVL
633	WT = UST(I) / PHIH
634	ENDIF
635	ZFUNC = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
636	EDYZ = KARMAN * WT * ZFUNC
637	! EDYZ = AMAX1(EDYZ,EDYZ0)
638	ENDIF
639	! PRINT *,I,K,TT(I,K)
640	!--------------------------------------------------------------------------
641	!
642	! RC = 0.257 * DZF ** 0.175
643	SS = ((US(I,K+1) - US(I,K)) 2 + (VS(I,K+1) - VS(I,K)) 2) &
644	/ (DZF * DZF) + 1.0E-9
645	GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
646	RI = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
647	! PRINT *,I,K,TT(I,K), RI, IMVDIF
648	!
649	!-- Adjustment to vert diff in Moist air
650	if(imvdif.eq.1)then
651	IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+ &
652	QIS(I,K+1)) .GT. 0.01E-3) THEN
653	QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
654	TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
655	XLV = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
656	ALPH = XLV * QMEAN / RD / TMEAN
657	CHI = XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
658	RI = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) * &
659	((CHI - ALPH) / (1.0 + CHI)))
660	ENDIF
661	endif
662
663	ZK = 0.4 * ZF(I,K)
664	SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
665
666	IF (RI .GE. RC) THEN
667	EDDYZ(I,K) = kzo !EDYZ0
668	ELSE IF (RI .GE. 0.0) THEN
669	EDDYZ(I,K) = kzo + SQRT(SS) * (1.- RI / RC) ** 2 * SQL
670	ELSE
671	EDDYZ(I,K) = kzo + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
672	ENDIF
673
674	IF(EDYZ.GT.EDDYZ(I,K)) THEN
675	EDDYZ(I,K) = EDYZ
676	ENDIF
677
678	EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
679	EDDYZ(I,K) = MAX(kzo,EDDYZ(I,K))
680
681	DENSF = 0.5 * (DENSX(I,K+1) + DENSX(I,K))
682
683	EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 * &
684	DTPBL * DSIGFI(K)*1E-6
685	ENDDO ! for I loop
686	ENDDO ! for k loop
687	!
688	DO I = its,ILX
689	EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
690	ENDDO
691	END SUBROUTINE EDDYX
692
693	!-------------------------------------------------------------------
694	SUBROUTINE ACM (DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
695	KLPBL, PBL, PBLSIG, MOL, UST, &
696	TST, QST, USTM, EDDYZ, DENSX, &
697	US, VS, THETA, QVS, QCS, QIS, &
698	UX, VX, THETAX, QVX, QCX, QIX, &
699	ids,ide, jds,jde, kds,kde, &
700	ims,ime, jms,jme, kms,kme, &
701	its,ite, jts,jte, kts,kte)
702	!**********************************************************************
703	! PBL model called the Asymmetric Convective Model, Version 2 (ACM2)
704	! -- See top of module for summary and references
705	!
706	!---- REVISION HISTORY:
707	! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
708	! JP and RG 8/2006 - updates
709	!
710	!**********************************************************************
711	! ARGUMENTS:
712	!-- DTPBL PBL time step
713	!-- PSTAR Psurf - Ptop in cb
714	!-- NOCONV If free convection =0, no; =1, yes
715	!-- SIGMAF Sigma for full layer
716	!-- DSIGH Sigma thickness
717	!-- DSIGHI Inverse of sigma thickness
718	!-- JX N-S index
719	!-- KLPBL PBL level at K index
720	!-- PBL PBL height in m
721	!-- PBLSIG Sigma level for PBL
722	!-- MOL Monin-Obukhov length in 1D form
723	!-- UST U* in 1D form
724	!-- TST Theta* in 1D form
725	!-- QST Q* in 1D form
726	!-- USTM U* for computation of momemtum flux
727	!-- EDDYZ eddy diffusivity KZ
728	!-- DENSX dry air density (kg/m^3)
729	!-- US U wind
730	!-- VS V wind
731	!-- THETA potential temperature
732	!-- QVS water vapor mixing ratio (Kg/Kg)
733	!-- QCS cloud mixing ratio (Kg/Kg)
734	!-- QIS ice mixing ratio (Kg/Kg)
735	!-- UX new U wind
736	!-- VX new V wind
737	!-- THETAX new potential temperature
738	!-- QVX new water vapor mixing ratio (Kg/Kg)
739	!-- QCX new cloud mixing ratio (Kg/Kg)
740	!-- QIX new ice mixing ratio (Kg/Kg)
741	!-----------------------------------------------------------------------
742
743	IMPLICIT NONE
744
745	!.......Arguments
746
747	!... Integer
748	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
749	ims,ime, jms,jme, kms,kme, &
750	its,ite, jts,jte, kts,kte, JX
751	INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV
752	INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL
753
754	!... Real
755	REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
756	REAL , INTENT(IN) :: DTPBL
757	REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, PBLSIG, &
758	MOL, TST, &
759	QST, USTM
760	REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGHI, DSIGH
761	REAL , DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
762	REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZ
763	REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, &
764	QVS, QCS, QIS, DENSX
765	REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX, &
766	QVX, QCX, QIX
767
768	!.......Local variables
769
770	!... Parameters
771	INTEGER, PARAMETER :: NSP = 6
772	!
773	!......ACM2 Parameters
774	! INTEGER, PARAMETER :: IFACM = 0
775	!
776	REAL, PARAMETER :: G1000 = 9.8 * 1.0E-3
777	REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP
778	REAL, PARAMETER :: KARMAN = 0.4
779
780	!... Integer
781	INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
782	INTEGER :: KCBLMX
783	INTEGER, DIMENSION( its:ite ) :: KCBL
784
785	!... Real
786	REAL :: G1000I, MBMAX, HOVL, MEDDY, MBAR
787	REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
788	REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, RAH, DTLIM
789	REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
790	REAL, DIMENSION( 1:NSP, its:ite ) :: FS, BCBOTN
791	REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS
792	REAL DELC
793	REAL, DIMENSION( 1:NSP,its:ite,kts:kte ) :: VCI
794
795	REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI !, Y
796	REAL, DIMENSION( 1:NSP,kts:kte ) :: DI, UI
797	!
798	!--Start Exicutable ----
799
800	ILX = ite
801	KL = kte
802	KLM = kte - 1
803
804	G1000I = 1.0 / G1000
805	KCBLMX = 0
806	MBMAX = 0.0
807
808	!---COMPUTE ACM MIXING RATE
809	DO I = its, ILX
810	DTLIM(I) = DTPBL
811	PSTARI(I) = 1.0 / PSTAR(I)
812	KCBL(I) = 1
813	FSACM(I) = 0.0
814
815	IF (NOCONV(I) .EQ. 1) THEN
816	KCBL(I) = KLPBL(I)
817	!-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
818	!--New couple ACM & EDDY-------------------------------------------------------------
819	HOVL = -PBL(I) / MOL(I)
820	FSACM(I) = 1./(1.+((KARMAN/(HOVL))*0.3333)/(0.72KARMAN))
821	MEDDY = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
822	MBAR = MEDDY * FSACM(I)
823	DO K = kts,KCBL(I)-1
824	EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
825	ENDDO
826	! if(i.eq.100.and.jx.eq.43) PRINT *,' Edy,meddy,mbar=', EDDYZ(I,1),MEDDY,MBAR
827	MBMAX = AMAX1(MBMAX,MBAR)
828	DO K = kts+1,KCBL(I)
829	MBARKS(K,I) = MBAR
830	MDWN(K,I) = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
831	ENDDO
832	MBARKS(1,I) = MBAR
833	MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
834	MDWN(KCBL(I)+1,I) = 0.0
835	ENDIF
836	ENDDO ! end of I loop
837
838	! if(NOCONV(20).EQ.1) print *,' KCBL,PBLSIG=',KCBL(20),PBLSIG(20),SIGMAF(KCBL(20)-1),DSIGHI(KCBL(20))
839	! do k = kts,klm
840	! if(NOCONV(20).EQ.1.and.MBAR.lt.1.0e-3)print *,' k,MBARKS,MDWN=',k,MBARKS(k,20),MDWN(k,20)
841	! enddo
842
843	DO K = kts,KLM
844	DO I = its,ILX
845	EKZ = EDDYZ(I,K) / DTPBL * DSIGHI(K)
846	DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
847	ENDDO
848	ENDDO
849
850	DO I = its,ILX
851	IF (NOCONV(I) .EQ. 1) THEN
852	KCBLMX = AMAX0(KLPBL(I),KCBLMX)
853	RZ = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
854	DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
855	ENDIF
856	ENDDO
857
858	DO I = its,ILX
859	ENDDO
860	!
861	DO K = kts,KL
862	DO I = its,ILX
863	VCI(1,I,K) = THETA(I,K)
864	VCI(2,I,K) = QVS(I,K)
865	VCI(3,I,K) = US(I,K)
866	VCI(4,I,K) = VS(I,K)
867	! -- Also mix cloud water and ice if necessary
868	! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN !!! Check other PBL models
869	VCI(5,I,K) = QCS(I,K)
870	VCI(6,I,K) = QIS(I,K)
871	ENDDO
872	ENDDO
873
874	NSPX=6
875
876	DO I = its,ILX
877	! RAH(I) = RA(I,J) + 3.976 / UST(I)
878	FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
879	FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
880	FM = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
881	WSPD = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
882	FS(3,I) = FM * US(I,1) / WSPD
883	FS(4,I) = FM * VS(I,1) / WSPD
884	FS(5,I) = 0.0
885	FS(6,I) = 0.0 ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
886	ENDDO
887	!
888	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
889	DO I = its,ILX
890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
891	NLP = INT(DTPBL / DTLIM(I) + 1.0)
892	DTS = (DTPBL / NLP)
893	DTRAT = DTS / DTPBL
894	DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP
895	!
896	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
897	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
898
899	!-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES
900
901	DO K = kts,KL
902	AI(K) = 0.0
903	BI(K) = 0.0
904	CI(K) = 0.0
905	EI(K) = 0.0
906	ENDDO
907
908	DO K = 2, KCBL(I)
909	EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
910	BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS
911	AI(K) = -CRANKP * MBARKS(K,I) * DTS
912	ENDDO
913
914	EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
915	AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT
916
917	DO K = KCBL(I)+1, KL
918	BI(K) = 1.0
919	ENDDO
920
921	DO K = 2,KL
922	XPLUS(K) = EDDYZ(I,K) * DSIGHI(K) * DTRAT
923	XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
924	CI(K) = - XMINUS(K) * CRANKP
925	EI(K) = EI(K) - XPLUS(K) * CRANKP
926	BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
927	ENDDO
928
929	IF (NOCONV(I) .EQ. 1) THEN
930	BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) * &
931	DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
932	ELSE
933	BI(1) = 1.0 + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
934	ENDIF
935
936
937	DO K = 1,KL
938	DO L = 1,NSPX
939	DI(L,K) = 0.0
940	ENDDO
941	ENDDO
942	!
943	!** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
944	DO K = 2,KCBL(I)
945	DO L = 1,NSPX
946	DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * &
947	VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) * &
948	MDWN(K+1,I) * VCI(L,I,K+1))
949	DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC
950	ENDDO
951	ENDDO
952
953	DO K = KCBL(I)+1, KL
954	DO L = 1,NSPX
955	DI(L,K) = VCI(L,I,K)
956	ENDDO
957	ENDDO
958
959	DO K = 2,KL
960	IF (K .EQ. KL) THEN
961	DO L = 1,NSPX
962	DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * &
963	(VCI(L,I,K) - VCI(L,I,K-1))
964	ENDDO
965	ELSE
966	DO L = 1,NSPX
967	DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * &
968	(VCI(L,I,K+1) - VCI(L,I,K)) - &
969	(1.0 - CRANKP) * XMINUS(K) * &
970	(VCI(L,I,K) - VCI(L,I,K-1))
971	ENDDO
972	ENDIF
973	ENDDO
974
975	IF (NOCONV(I) .EQ. 1) THEN
976	DO L = 1,NSPX
977	F1 = -G1000I * (MBARKS(1,I) * &
978	(PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) - &
979	MDWN(2,I) * VCI(L,I,2) * DSIGH(2))
980
981	DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP) &
982	* F1) * DSIGHI(1) * DTS
983	ENDDO
984	ELSE
985	DO L = 1,NSPX
986	DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
987	ENDDO
988	ENDIF
989	DO L = 1,NSPX
990	DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1) &
991	* DTRAT * (VCI(L,I,2) - VCI(L,I,1))
992	ENDDO
993	IF ( NOCONV(I) .EQ. 1 ) THEN
994	CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
995	ELSE
996	CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
997	END IF
998	!
999	!-- COMPUTE NEW THETAV AND Q
1000	DO K = 1,KL
1001	DO L = 1,NSPX
1002	VCI(L,I,K) = UI(L,K)
1003	ENDDO
1004	ENDDO
1005	!
1006	ENDDO ! END I LOOP
1007	ENDDO ! END SUB TIME LOOP
1008	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1009	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1010
1011	!
1012	DO K = kts,KL
1013	DO I = its,ILX
1014	THETAX(I,K) = VCI(1,I,K)
1015	QVX(I,K) = VCI(2,I,K)
1016	UX(I,K) = VCI(3,I,K)
1017	VX(I,K) = VCI(4,I,K)
1018	! IF( (I.EQ.61 .OR. I.EQ.15) .AND. JX.EQ.22) THEN
1019	! PRINT *,'TEST -->',I,K,THETAX(I,K)
1020	! ENDIF
1021	ENDDO
1022	ENDDO
1023
1024	DO K = kts,KL
1025	DO I = its,ILX
1026	QCX(I,K) = VCI(5,I,K)
1027	QIX(I,K) = VCI(6,I,K)
1028	ENDDO
1029	ENDDO
1030	! ENDIF
1031
1032	END SUBROUTINE ACM
1033	!--------------------------------------------------------
1034
1035	SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
1036
1037	IMPLICIT NONE
1038	!
1039	!-- Bordered band diagonal matrix solver for ACM2
1040
1041	!-- ACM2 Matrix is in this form:
1042	! B1 E1
1043	! A2 B2 E2
1044	! A3 C3 B3 E3
1045	! A4 C4 B4 E4
1046	! A5 C5 B5 E5
1047	! A6 C6 B6
1048
1049	!--Upper Matrix is
1050	! U11 U12
1051	! U22 U23
1052	! U33 U34
1053	! U44 U45
1054	! U55 U56
1055	! U66
1056
1057	!--Lower Matrix is:
1058	! 1
1059	! L21 1
1060	! L31 L32 1
1061	! L41 L42 L43 1
1062	! L51 L52 L53 L54 1
1063	! L61 L62 L63 L64 L65 1
1064	!---------------------------------------------------------
1065	!...Arguments
1066	INTEGER, INTENT(IN) :: KL
1067	INTEGER, INTENT(IN) :: NSP
1068	REAL A(KL),B(KL),E(KL)
1069	REAL C(KL),D(NSP,KL),X(NSP,KL)
1070
1071	!...Locals
1072	REAL Y(NSP,KL),AIJ,SUM
1073	REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
1074	INTEGER I,J,V
1075
1076	!-- Define Upper and Lower matrices
1077	L(1,1) = 1.
1078	! U(1,1) = B(1)
1079	UII(1) = B(1)
1080	RUII(1) = 1./UII(1)
1081	DO I = 2, KL
1082	L(I,I) = 1.
1083	L(I,1) = A(I)/B(1)
1084	! U(I-1,I) =E(I-1)
1085	UIIP1(I-1)=E(I-1)
1086	IF(I.GE.3) THEN
1087	DO J = 2,I-1
1088	IF(I.EQ.J+1) THEN
1089	AIJ = C(I)
1090	ELSE
1091	AIJ = 0.
1092	ENDIF
1093	L(I,J) = (AIJ-L(I,J-1)*E(J-1))/ &
1094	(B(J)-L(J,J-1)*E(J-1))
1095	ENDDO
1096	ENDIF
1097	ENDDO
1098
1099	DO I = 2,KL
1100	! U(I,I) = B(I)-L(I,I-1)*E(I-1)
1101	UII(I) = B(I)-L(I,I-1)*E(I-1)
1102	RUII(I) = 1./UII(I)
1103	ENDDO
1104
1105	!-- Forward sub for Ly=d
1106	DO V= 1, NSP
1107	Y(V,1) = D(V,1)
1108	DO I=2,KL
1109	SUM = D(V,I)
1110	DO J=1,I-1
1111	SUM = SUM - L(I,J)*Y(V,J)
1112	ENDDO
1113	Y(V,I) = SUM
1114	ENDDO
1115	ENDDO
1116
1117	!-- Back sub for Ux=y
1118
1119	DO V= 1, NSP
1120	! X(V,KL) = Y(V,KL)/U(KL,KL)
1121	X(V,KL) = Y(V,KL)*RUII(KL)
1122	ENDDO
1123	DO I = KL-1,1,-1
1124	DO V= 1, NSP
1125	! X(V,I) = (Y(V,I)-U(I,I+1)*X(V,I+1))/U(I,I)
1126	X(V,I) = (Y(V,I)-UIIP1(I)X(V,I+1))RUII(I)
1127	ENDDO
1128	ENDDO
1129
1130	END SUBROUTINE MATRIX
1131	!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1132	SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
1133	!-----------------------------------------------------------------------
1134
1135	! FUNCTION:
1136	! Solves tridiagonal system by Thomas algorithm.
1137	! The associated tri-diagonal system is stored in 3 arrays
1138	! D : diagonal
1139	! L : sub-diagonal
1140	! U : super-diagonal
1141	! B : right hand side function
1142	! X : return solution from tridiagonal solver
1143
1144	! [ D(1) U(1) 0 0 0 ... 0 ]
1145	! [ L(2) D(2) U(2) 0 0 ... . ]
1146	! [ 0 L(3) D(3) U(3) 0 ... . ]
1147	! [ . . . . . ] X(i) = B(i)
1148	! [ . . . . 0 ]
1149	! [ . . . . ]
1150	! [ 0 L(n) D(n) ]
1151
1152	!-----------------------------------------------------------------------
1153
1154	IMPLICIT NONE
1155
1156	! Arguments:
1157
1158	INTEGER, INTENT(IN) :: KL
1159	INTEGER, INTENT(IN) :: NSP
1160
1161	REAL L( KL ) ! subdiagonal
1162	REAL D(KL) ! diagonal
1163	REAL U( KL ) ! superdiagonal
1164	REAL B(NSP,KL ) ! R.H. side
1165	REAL X( NSP,KL ) ! solution
1166
1167	! Local Variables:
1168
1169	REAL GAM( KL )
1170	REAL BET
1171	INTEGER V, K
1172
1173	! Decomposition and forward substitution:
1174	BET = 1.0 / D( 1 )
1175	DO V = 1, NSP
1176	X( V,1 ) = BET * B(V,1 )
1177	ENDDO
1178
1179	DO K = 2, KL
1180	GAM(K ) = BET * U( K-1 )
1181	BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
1182	DO V = 1, NSP
1183	X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
1184	ENDDO
1185	END DO
1186
1187	! Back-substitution:
1188
1189	DO K = KL - 1, 1, -1
1190	DO V = 1, NSP
1191	X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
1192	END DO
1193	ENDDO
1194
1195	END SUBROUTINE TRI
1196
1197	END MODULE module_bl_acm
1198

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format