Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_acm.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

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

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

Download in other formats:

Original Format