Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_mrf.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:

File size: 80.1 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE module_bl_mrf
4
5	CONTAINS
6
7	!-------------------------------------------------------------------
8	SUBROUTINE MRF(U3D,V3D,TH3D,T3D,QV3D,QC3D,P3D,PI3D, &
9	RUBLTEN,RVBLTEN,RTHBLTEN, &
10	RQVBLTEN,RQCBLTEN, &
11	CP,G,ROVCP,R,ROVG, &
12	dz8w,z,XLV,RV,PSFC, &
13	p1000mb, &
14	ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
15	XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
16	DT,DTMIN,KPBL2D, &
17	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
18	flag_qi, &
19	ids,ide, jds,jde, kds,kde, &
20	ims,ime, jms,jme, kms,kme, &
21	its,ite, jts,jte, kts,kte, &
22	! Optional
23	QI3D,RQIBLTEN, &
24	regime )
25	!-------------------------------------------------------------------
26	IMPLICIT NONE
27	!-------------------------------------------------------------------
28	!-- U3D 3D u-velocity interpolated to theta points (m/s)
29	!-- V3D 3D v-velocity interpolated to theta points (m/s)
30	!-- TH3D 3D potential temperature (K)
31	!-- T3D temperature (K)
32	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
33	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
34	!-- QI3D 3D ice mixing ratio (Kg/Kg)
35	!-- P3D 3D pressure (Pa)
36	!-- PI3D 3D exner function (dimensionless)
37	!-- rr3D 3D dry air density (kg/m^3)
38	!-- RUBLTEN U tendency due to
39	! PBL parameterization (m/s^2)
40	!-- RVBLTEN V tendency due to
41	! PBL parameterization (m/s^2)
42	!-- RTHBLTEN Theta tendency due to
43	! PBL parameterization (K/s)
44	!-- RQVBLTEN Qv tendency due to
45	! PBL parameterization (kg/kg/s)
46	!-- RQCBLTEN Qc tendency due to
47	! PBL parameterization (kg/kg/s)
48	!-- RQIBLTEN Qi tendency due to
49	! PBL parameterization (kg/kg/s)
50	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
51	!-- G acceleration due to gravity (m/s^2)
52	!-- ROVCP R/CP
53	!-- R gas constant for dry air (J/kg/K)
54	!-- ROVG R/G
55	!-- dz8w dz between full levels (m)
56	!-- z height above sea level (m)
57	!-- XLV latent heat of vaporization (J/kg)
58	!-- RV gas constant for water vapor (J/kg/K)
59	!-- PSFC pressure at the surface (Pa)
60	!-- ZNT roughness length (m)
61	!-- UST u* in similarity theory (m/s)
62	!-- ZOL z/L height over Monin-Obukhov length
63	!-- HOL PBL height over Monin-Obukhov length
64	!-- PBL PBL height (m)
65	!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
66	!-- PSIM similarity stability function for momentum
67	!-- PSIH similarity stability function for heat
68	!-- XLAND land mask (1 for land, 2 for water)
69	!-- HFX upward heat flux at the surface (W/m^2)
70	!-- QFX upward moisture flux at the surface (kg/m^2/s)
71	!-- TSK surface temperature (K)
72	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
73	!-- WSPD wind speed at lowest model level (m/s)
74	!-- BR bulk Richardson number in surface layer
75	!-- DT time step (s)
76	!-- DTMIN time step (minute)
77	!-- rvovrd R_v divided by R_d (dimensionless)
78	!-- SVP1 constant for saturation vapor pressure (kPa)
79	!-- SVP2 constant for saturation vapor pressure (dimensionless)
80	!-- SVP3 constant for saturation vapor pressure (K)
81	!-- SVPT0 constant for saturation vapor pressure (K)
82	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
83	!-- EP2 constant for specific humidity calculation
84	!-- KARMAN Von Karman constant
85	!-- EOMEG angular velocity of earth's rotation (rad/s)
86	!-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
87	!-- ids start index for i in domain
88	!-- ide end index for i in domain
89	!-- jds start index for j in domain
90	!-- jde end index for j in domain
91	!-- kds start index for k in domain
92	!-- kde end index for k in domain
93	!-- ims start index for i in memory
94	!-- ime end index for i in memory
95	!-- jms start index for j in memory
96	!-- jme end index for j in memory
97	!-- kms start index for k in memory
98	!-- kme end index for k in memory
99	!-- its start index for i in tile
100	!-- ite end index for i in tile
101	!-- jts start index for j in tile
102	!-- jte end index for j in tile
103	!-- kts start index for k in tile
104	!-- kte end index for k in tile
105	!-------------------------------------------------------------------
106
107	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
108	ims,ime, jms,jme, kms,kme, &
109	its,ite, jts,jte, kts,kte
110
111	!
112	REAL, INTENT(IN ) :: P1000mb
113	REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
114	ROVG,R,XLV,RV
115
116	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
117	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
118
119	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
120	INTENT(IN ) :: QV3D, &
121	QC3D, &
122	P3D, &
123	PI3D, &
124	TH3D, &
125	T3D, &
126	dz8w, &
127	z
128	!
129	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
130	INTENT(INOUT) :: RUBLTEN, &
131	RVBLTEN, &
132	RTHBLTEN, &
133	RQVBLTEN, &
134	RQCBLTEN
135
136	REAL, DIMENSION( ims:ime, jms:jme ) , &
137	INTENT(IN ) :: XLAND, &
138	HFX, &
139	QFX
140
141	REAL, DIMENSION( ims:ime, jms:jme ) , &
142	INTENT(INOUT) :: HOL, &
143	UST, &
144	PBL, &
145	ZNT
146
147	LOGICAL, INTENT(IN) :: FLAG_QI
148	!
149	!m The following 5 variables are changed to memory size from tile size--
150	!
151	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
152	PSIM, &
153	PSIH
154
155	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
156	WSPD
157
158	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
159	GZ1OZ0, &
160	BR
161
162	REAL, DIMENSION( ims:ime, jms:jme ) , &
163	INTENT(IN ) :: PSFC
164
165	REAL, DIMENSION( ims:ime, jms:jme ) , &
166	INTENT(IN ) :: TSK
167
168	REAL, DIMENSION( ims:ime, jms:jme ) , &
169	INTENT(INOUT) :: ZOL
170
171	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
172	INTENT(OUT ) :: KPBL2D
173
174	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
175	INTENT(IN ) :: U3D, &
176	V3D
177	!
178	! Optional
179	!
180	REAL, DIMENSION( ims:ime, jms:jme ) , &
181	OPTIONAL , &
182	INTENT(INOUT) :: REGIME
183
184	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
185	INTENT(INOUT) :: RQIBLTEN
186
187	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
188	OPTIONAL , &
189	INTENT(IN ) :: QI3D
190
191	! LOCAL VARS
192	REAL, DIMENSION( its:ite, kts:kte ) :: dz8w2d, &
193	z2d
194
195
196	INTEGER :: I,J,K,NK
197
198	!
199	DO J=jts,jte
200	DO k=kts,kte
201	NK=kme-k
202	DO i=its,ite
203	dz8w2d(I,K) = dz8w(i,NK,j)
204	z2d(I,K) = z(i,NK,j)
205	ENDDO
206	ENDDO
207
208
209	CALL MRF2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j), &
210	QV3D(ims,kms,j),QC3D(ims,kms,j), &
211	P3D(ims,kms,j),RUBLTEN(ims,kms,j),RVBLTEN(ims,kms,j),&
212	RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j), &
213	RQCBLTEN(ims,kms,j), &
214	p1000mb, &
215	CP,G,ROVCP,R,ROVG, &
216	dz8w2d,z2d,XLV,Rv, &
217	PSFC(ims,j),ZNT(ims,j), &
218	UST(ims,j),ZOL(ims,j), &
219	HOL(ims,j),PBL(ims,j),PSIM(ims,j), &
220	PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j), &
221	TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j), &
222	DT,DTMIN,KPBL2D(ims,j), &
223	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
224	flag_qi, &
225	ids,ide, jds,jde, kds,kde, &
226	ims,ime, jms,jme, kms,kme, &
227	its,ite, jts,jte, kts,kte, &
228	!optional
229	QI2DTEN=RQIBLTEN(ims,kms,j), &
230	REGIME=REGIME(ims,j),QI2D=QI3D(ims,kms,j) )
231
232
233	DO k=kts,kte
234	DO i=its,ite
235	RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
236	ENDDO
237	ENDDO
238
239	ENDDO
240
241	END SUBROUTINE MRF
242
243	!-------------------------------------------------------------------
244	SUBROUTINE MRF2D(J,U2D,V2D,T2D,QV2D,QC2D, P2D, &
245	U2DTEN,V2DTEN,T2DTEN, &
246	QV2DTEN,QC2DTEN, &
247	p1000mb, &
248	CP,G,ROVCP,R,ROVG, &
249	dz8w2d,z2d,XLV,RV,PSFCPA, &
250	ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
251	XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
252	DT,DTMIN,KPBL1D, &
253	SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
254	flag_qi, &
255	ids,ide, jds,jde, kds,kde, &
256	ims,ime, jms,jme, kms,kme, &
257	its,ite, jts,jte, kts,kte, &
258	! optional
259	regime, qi2d, QI2DTEN )
260	!-------------------------------------------------------------------
261	IMPLICIT NONE
262	!-------------------------------------------------------------------
263	! BASED ON THE "COUNTERGRADIENT" TRANSPORT TERM OF TROEN
264	! AND MAHRT (1986) FOR THE UNSTABLE PBL.
265	! THIS ROUTINE USES AN IMPLICIT APPROACH FOR VERTICAL FLUX
266	! DIVERGENCE AND DOES NOT REQUIRE "MITER" TIMESTEPS.
267	! IT INCLUDES VERTICAL DIFFUSION IN THE STABLE ATMOSPHERE
268	! AND MOIST VERTICAL DIFFUSION IN CLOUDS.
269	! SURFACE FLUXES CALCULATED AS IN HIRPBL.
270	! 5-LAYER SOIL MODEL OPTION REQUIRED IN SLAB DUE TO LONG TIMESTEP
271	!
272	! CODED BY SONG-YOU HONG (NCEP), IMPLEMENTED BY JIMY DUDHIA (NCAR)
273	! FALL 1996
274	!
275	! REFERENCES:
276	!
277	! HONG AND PAN (1996), MON. WEA. REV.
278	! TROEN AND MAHRT (1986), BOUNDARY LAYER MET.
279	!
280	! CHANGES:
281	! INCREASE RLAM FROM 30 TO 150, AND CHANGE FREE ATMOSPHERE
282	! STABILITY FUNCTION TO INCREASE VERTICAL DIFFUSION
283	! (HONG, JUNE 1997)
284	!
285	! PUT LOWER LIMIT ON PSI FOR STABLE CONDITIONS. THIS WILL
286	! PREVENT FLUXES FROM BECOMING TOO SMALL (DUDHIA, OCTOBER 1997)
287	!
288	! CORRECTION TO REGIME CALCULATION. THIS WILL ALLOW POINTS IN
289	! REGIME 4 MUCH MORE FREQUENTLY GIVING LARGER SURFACE FLUXES
290	! REGIME 3 NO LONGER USES HOL < 1.5 OR THVX LAPSE-RATE CHECK
291	! IN MRF SCHEME. THIS WILL MAKE REGIME 3 MUCH LESS FREQUENT.
292	!
293	! ADD SURFACE PRESSURE, PS(I), ARRAY FOR EFFICIENCY
294	!
295	! FIX FOR PROBLEM WITH THIN LAYERS AND HIGH ROUGHNESS
296	!
297	! CHARNOCK CONSTANT NOW COMES FROM NAMELIST (DEFAULT SAME)
298	!
299	!-------------------------------------------------------------------
300
301	REAL RLAM,PRMIN,PRMAX,XKZMIN,XKZMAX,RIMIN,BRCR, &
302	CFAC,PFAC,SFCFRAC,CKZ,ZFMIN,APHI5,APHI16,GAMCRT, &
303	GAMCRQ,XKA,PRT
304
305	PARAMETER (RLAM=150.,PRMIN=0.5,PRMAX=4.)
306	PARAMETER (XKZMIN=0.01,XKZMAX=1000.,RIMIN=-100.)
307	PARAMETER (BRCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
308	PARAMETER (CKZ=0.001,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
309	PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)
310	PARAMETER (XKA=2.4E-5)
311	PARAMETER (PRT=1.)
312	!
313	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
314	ims,ime, jms,jme, kms,kme, &
315	its,ite, jts,jte, kts,kte, &
316	J
317	!
318	LOGICAL, INTENT(IN) :: FLAG_QI
319	!
320	REAL, INTENT(IN ) :: P1000mb
321	REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
322	ROVG,R,XLV,RV
323
324	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
325	REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
326
327	REAL, DIMENSION( ims:ime, kms:kme ) , &
328	INTENT(IN ) :: QV2D, &
329	QC2D, &
330	P2D, &
331	T2D
332	!
333	REAL, DIMENSION( ims:ime, kms:kme ) , &
334	INTENT(INOUT) :: U2DTEN, &
335	V2DTEN, &
336	T2DTEN, &
337	QV2DTEN, &
338	QC2DTEN
339
340	REAL, DIMENSION( ims:ime ) , &
341	INTENT(INOUT) :: HOL, &
342	UST, &
343	PBL, &
344	ZNT
345
346	REAL, DIMENSION( ims:ime ) , &
347	INTENT(IN ) :: XLAND, &
348	HFX, &
349	QFX
350	!
351	!m The following 5 are changed to memory size---
352	!
353	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: PSIM, &
354	PSIH
355
356	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: WSPD
357
358	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: GZ1OZ0, &
359	BR
360
361	REAL, DIMENSION( ims:ime ) , &
362	INTENT(IN ) :: PSFCPA
363
364	REAL, DIMENSION( ims:ime ) , &
365	INTENT(IN ) :: TSK
366
367	REAL, DIMENSION( ims:ime ) , &
368	INTENT(INOUT) :: ZOL
369
370	INTEGER, DIMENSION( ims:ime ) , &
371	INTENT(OUT ) :: KPBL1D
372
373	REAL, DIMENSION( ims:ime, kms:kme ) , &
374	INTENT(IN ) :: U2D, &
375	V2D
376
377	! MODULE-LOCAL VARIABLES (DEFINED IN SUBROUTINE MRF)
378	!
379	REAL, DIMENSION( its:ite, kts:kte ) , &
380	INTENT(IN) :: dz8w2d, &
381	z2d
382	!
383	!
384	! Optional
385	!
386	REAL, DIMENSION( ims:ime ) , &
387	OPTIONAL , &
388	INTENT(INOUT) :: REGIME
389
390	REAL, DIMENSION( ims:ime, kms:kme ) , &
391	OPTIONAL , &
392	INTENT(IN ) :: QI2D
393
394	REAL, DIMENSION( ims:ime, kms:kme ) , &
395	OPTIONAL , &
396	INTENT(INOUT) :: QI2DTEN
397
398	! LOCAL VARS
399
400	REAL, DIMENSION( its:ite, kts:kte+1 ) :: ZQ
401
402	REAL, DIMENSION( its:ite, kts:kte ) :: &
403	UX,VX,QX, &
404	QCX,THX,THVX, &
405	DZQ,DZA, &
406	TTNP,QTNP, &
407	QCTNP,ZA, &
408	UXS,VXS, &
409	THXS,QXS, &
410	QCXS,QIX, &
411	QITNP,QIXS, &
412	UTNP,VTNP
413	!
414	REAL, DIMENSION( its:ite ) :: QIXSV,RHOX, &
415	WSPD1,GOVRTH, &
416	PBL0,THXSV, &
417	UXSV,VXSV, &
418	QXSV,QCXSV, &
419	QGH,TGDSA,PS
420
421	INTEGER :: ILXM,JLXM,KL, &
422	KLM,KLP1,KLPBL
423	!
424	INTEGER, DIMENSION( its:ite ) :: KPBL,KPBL0
425	!
426	REAL, DIMENSION( its:ite, kts:kte ) :: SCR3,SCR4
427	!
428	REAL, DIMENSION( its:ite ) :: DUM1, &
429	XKZMKL
430	!
431	REAL, DIMENSION( its:ite ) :: ZL1,THERMAL, &
432	WSCALE,HGAMT, &
433	HGAMQ,BRDN, &
434	BRUP,PHIM, &
435	PHIH,CPM, &
436	DUSFC,DVSFC, &
437	DTSFC,DQSFC
438
439	!
440	REAL, DIMENSION( its:ite, kts:kte ) :: XKZM,XKZH, &
441	A1,A2, &
442	AD,AU, &
443	TX
444	!
445	REAL, DIMENSION( its:ite, kts:kte ) :: AL
446	!
447	LOGICAL, DIMENSION( its:ite ) :: PBLFLG, &
448	SFCFLG, &
449	STABLE
450	!
451	REAL, DIMENSION( its:ite ) :: THGB
452
453	INTEGER :: N,I,K,KK,L,NZOL,IMVDIF
454
455	INTEGER :: JBGN,JEND,IBGN,IEND,NK
456
457	REAL :: ZOLN,X,Y,CONT,CONQ,CONW,PL,THCON,TVCON,E1,DTSTEP
458	REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL
459	REAL :: DTTHX,PSIX,DTG,PSIQ,USTM
460	REAL :: DT4,RDT,SPDK2,FM,FH,HOL1,GAMFAC,VPERT,PRNUM
461	REAL :: ZFAC,XKZO,SS,RI,QMEAN,TMEAN,ALPH,CHI,ZK,RL2,DK,SRI
462	REAL :: BRINT,DTODSD,DSIG,RDZ,DSDZT,DSDZQ,DSDZ2,TTEND,QTEND
463	REAL :: UTEND,VTEND,QCTEND,QITEND,TGC,DTODSU
464
465	!----------------------------------------------------------------------
466
467	KLPBL=1
468	KL=kte
469	ILXM=ite-1
470	JLXM=jte-1
471	KLM=kte-1
472	KLP1=kte+1
473	!
474	CONT=1000.*CP/G
475	CONQ=1000.*XLV/G
476	CONW=1000./G
477
478	!-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
479
480	IMVDIF=1
481
482	! DO i=its,ite
483	!!PS PSFC cmb
484	! PSFC(I)=PSFCPA(I)/1000.
485	! ENDDO
486
487
488	!
489	!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
490	!
491	DO 5 I=its,ite
492	TGDSA(I)=TSK(I)
493	! PS PSFC cmb
494	PS(I)=PSFCPA(I)/1000.
495	! THGB(I)=TSK(I)(100./PS(I))*ROVCP
496	THGB(I)=TSK(I)(P1000mb/PSFCPA(I))*ROVCP
497	5 CONTINUE
498	!
499	!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
500	! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
501	!
502	! * NOTE *
503	! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
504	! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
505	! TENDENCIES.
506	!
507	DO 24 K=kts,kte
508	NK=kme-K
509	DO 24 I=its,ite
510	UX(I,K)=U2D(I,NK)
511	VX(I,K)=V2D(I,NK)
512	24 CONTINUE
513	!
514	!.....SCR3(I,K) STORE TEMPERATURE,
515	! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
516	!
517	DO 30 K=kts,kte
518	NK=kme-K
519	DO 30 I=its,ite
520	! PL cmb
521	PL=P2D(I,NK)/1000.
522	SCR3(I,K)=T2D(I,NK)
523	! THCON=(100./PL)**ROVCP
524	THCON=(P1000mb/(PL1000.))*ROVCP
525	THX(I,K)=SCR3(I,K)*THCON
526	TX(I,K)=SCR3(I,K)
527	SCR4(I,K)=SCR3(I,K)
528	THVX(I,K)=THX(I,K)
529	QX(I,K)=0.
530	30 CONTINUE
531	!
532	DO I=its,ite
533	QGH(i)=0.
534	CPM(i)=CP
535	ENDDO
536	!
537	! IF(IDRY.EQ.1)GOTO 80
538	DO 50 K=kts,kte
539	NK=kme-K
540	DO 50 I=its,ite
541	QX(I,K)=QV2D(I,NK)
542	TVCON=(1.+EP1*QX(I,K))
543	THVX(I,K)=THX(I,K)*TVCON
544	SCR4(I,K)=SCR3(I,K)*TVCON
545	50 CONTINUE
546	!
547	DO 60 I=its,ite
548	E1=SVP1EXP(SVP2(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
549	QGH(I)=EP2*E1/(PS(I)-E1)
550	CPM(I)=CP(1.+0.8QX(I,KL))
551	60 CONTINUE
552	!
553	! IF(IMOIST.EQ.1)GOTO 80
554	DO 70 K=kts,kte
555	NK=kme-K
556	DO 70 I=its,ite
557	QCX(I,K)=QC2D(I,NK)
558	70 CONTINUE
559
560	IF (flag_QI .AND. PRESENT( QI2D ) ) THEN
561	DO K=kts,kte
562	NK=kme-K
563	DO I=its,ite
564	QIX(I,K)=QI2D(I,NK)
565	ENDDO
566	ENDDO
567	ELSE
568	DO K=kts,kte
569	NK=kme-K
570	DO I=its,ite
571	QIX(I,K)=0.
572	ENDDO
573	ENDDO
574	ENDIF
575
576	80 CONTINUE
577
578	!
579	!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
580	! LEVEL, AND THE LAYER THICKNESSES.
581	!
582	DO 90 I=its,ite
583	ZQ(I,KLP1)=0.
584	RHOX(I)=PS(I)1000./(RSCR4(I,KL))
585	90 CONTINUE
586	!
587	DO 110 KK=kts,kte
588	K=kme-KK
589	DO 100 I=its,ite
590	DUM1(I)=ZQ(I,K+1)
591	100 CONTINUE
592	!
593	DO 110 I=its,ite
594	ZQ(I,K)=dz8w2d(I,K)+DUM1(I)
595	110 CONTINUE
596	!
597	DO 120 K=kts,kte
598	DO 120 I=its,ite
599	ZA(I,K)=0.5*(ZQ(I,K)+ZQ(I,K+1))
600	DZQ(I,K)=ZQ(I,K)-ZQ(I,K+1)
601	120 CONTINUE
602	!
603	DO 130 K=kts,kte-1
604	DO 130 I=its,ite
605	DZA(I,K)=ZA(I,K)-ZA(I,K+1)
606	130 CONTINUE
607
608	DTSTEP=DT
609	!
610	DO 160 I=its,ite
611	GOVRTH(I)=G/THX(I,KL)
612	160 CONTINUE
613	!
614	!-----INITIALIZE VERTICAL TENDENCIES AND
615	!
616	DO I=its,ite
617	DO K=kts,kte
618	UTNP(i,k)=0.
619	VTNP(i,k)=0.
620	TTNP(i,k)=0.
621	ENDDO
622	ENDDO
623	!
624	! IF(IDRY.EQ.1)GOTO 250
625	DO 230 K=kts,kte
626	DO 230 I=its,ite
627	QTNP(I,K)=0.
628	230 CONTINUE
629	!
630	! IF(IMOIST.EQ.1)GOTO 250
631	DO 240 K=kts,kte
632	DO 240 I=its,ite
633	QCTNP(I,K)=0.
634	QITNP(I,K)=0.
635	240 CONTINUE
636
637	250 CONTINUE
638	!
639	!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
640	! AKB(1976), EQ(12).
641
642	! DO 260 I=its,ite
643	! GZ1OZ0(I)=ALOG(ZA(I,KL)/ZNT(I))
644	! IF((XLAND(I)-1.5).GE.0)THEN
645	! ZL=ZNT(I)
646	! ELSE
647	! ZL=0.01
648	! ENDIF
649	! WSPD(I)=SQRT(UX(I,KL)UX(I,KL)+VX(I,KL)VX(I,KL))
650	! TSKV=THGB(I)(1.+EP1QGH(I)*MAVAIL(I))
651	! DTHVDZ=(THVX(I,KL)-TSKV)
652	! IF(-DTHVDZ.GE.0)THEN
653	! DTHVM=-DTHVDZ
654	! ELSE
655	! DTHVM=0.
656	! ENDIF
657	! VCONV=VCONVC*SQRT(DTHVM)
658	! WSPD(I)=SQRT(WSPD(I)WSPD(I)+VCONVVCONV)
659	! WSPD(I)=AMAX1(WSPD(I),1.)
660	! BR(I)=GOVRTH(I)ZA(I,KL)DTHVDZ/(WSPD(I)*WSPD(I))
661	! 260 CONTINUE
662
663	!!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
664	!!
665	!! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
666	!! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
667	!!
668	!! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
669	!!
670	!! 1. BR .GE. 0.2;
671	!! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
672	!!
673	!! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
674	!! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
675	!! (REGIME=2),
676	!!
677	!! 3. BR .EQ. 0.0
678	!! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
679	!!
680	!! 4. BR .LT. 0.0
681	!! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
682	!!
683	!!-----
684	!
685	! DO 320 I=its,ite
686	!!----
687	!!-- REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
688	!!-- IF(BR(I).LT.0..AND.HOL(I).GT.1.5)GOTO 310
689	!
690	! IF(BR(I).LT.0.)GOTO 310
691	!!
692	!!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
693	!!
694	! IF(BR(I).LT.0.2)GOTO 270
695	! REGIME(I)=1.
696	! PSIM(I)=-10.*GZ1OZ0(I)
697	!! LOWER LIMIT ON PSI IN STABLE CONDITIONS
698	! PSIM(I)=AMAX1(PSIM(I),-10.)
699	! PSIH(I)=PSIM(I)
700	! HOL(I)=0.0
701	! PBL(I)=0.0
702	! GOTO 320
703	!!
704	!!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
705	!!
706	! 270 IF(BR(I).EQ.0.0)GOTO 280
707	! REGIME(I)=2.
708	! PSIM(I)=-5.0BR(I)GZ1OZ0(I)/(1.1-5.0*BR(I))
709	!! LOWER LIMIT ON PSI IN STABLE CONDITIONS
710	! PSIM(I)=AMAX1(PSIM(I),-10.)
711	!!.....AKB(1976), EQ(16).
712	! PSIH(I)=PSIM(I)
713	! HOL(I)=0.0
714	! PBL(I)=0.0
715	! GOTO 320
716	!!
717	!!-----CLASS 3; FORCED CONVECTION:
718	!!
719	! 280 REGIME(I)=3.
720	! PSIM(I)=0.0
721	! PSIH(I)=PSIM(I)
722	!
723	!! special use kte instead of kme
724	!
725	! DO 290 KK=kts,kte-1
726	! K=kte-KK
727	! IF(THVX(I,K).GT.THVX(I,KL))GOTO 300
728	! 290 CONTINUE
729	! STOP 290
730	! 300 PBL(I)=ZQ(I,K+1)
731	! IF(UST(I).LT.0.01)THEN
732	! ZOL(I)=BR(I)*GZ1OZ0(I)
733	! ELSE
734	! ZOL(I)=KARMANGOVRTH(I)ZA(I,KL)MOL(I,J)/(UST(I)UST(I))
735	! ENDIF
736	! HOL(I)=-ZOL(I)*PBL(I)/ZA(I,KL)
737	! GOTO 320
738	!
739	!!-----CLASS 4; FREE CONVECTION:
740	!
741	!! 310 IF(THVX(I,KLM).GT.THVX(I,KL))GOTO 280
742	!
743	! 310 CONTINUE
744	! REGIME(I)=4.
745	! IF(UST(I).LT.0.01)THEN
746	! ZOL(I)=BR(I)*GZ1OZ0(I)
747	! ELSE
748	! ZOL(I)=KARMANGOVRTH(I)ZA(I,KL)MOL(I,J)/(UST(I)UST(I))
749	! ENDIF
750	! ZOL(I)=AMIN1(ZOL(I),0.)
751	! ZOL(I)=AMAX1(ZOL(I),-9.9999)
752	! NZOL=INT(-ZOL(I)*100.)
753	! RZOL=-ZOL(I)*100.-NZOL
754	! PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
755	! PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
756	!!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
757	!!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
758	! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
759	! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
760	! 320 CONTINUE
761
762	!-----COMPUTE THE FRICTIONAL VELOCITY:
763	! ZA(1982) EQS(2.60),(2.61).
764
765	DO 330 I=its,ite
766	DTG=THX(I,KL)-THGB(I)
767	PSIX=GZ1OZ0(I)-PSIM(I)
768	IF((XLAND(I)-1.5).GE.0)THEN
769	ZL=ZNT(I)
770	ELSE
771	ZL=0.01
772	ENDIF
773	PSIQ=ALOG(KARMANUST(I)ZA(I,KL)/XKA+ZA(I,KL)/ZL)-PSIH(I)
774	UST(I)=KARMAN*WSPD(I)/PSIX
775	!
776	USTM=AMAX1(UST(I),0.1)
777	IF((XLAND(I)-1.5).GE.0)THEN
778	UST(I)=UST(I)
779	ELSE
780	UST(I)=USTM
781	ENDIF
782	! MOL(I,J)=KARMAN*DTG/(GZ1OZ0(I)-PSIH(I))/PRT
783	330 CONTINUE
784	!
785	DO 420 I=its,ite
786	WSPD1(I)=SQRT(UX(I,KL)UX(I,KL)+VX(I,KL)VX(I,KL))+1.E-9
787	420 CONTINUE
788	!
789	!---- COMPUTE VERTICAL DIFFUSION
790	!
791	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
792	! COMPUTE PRELIMINARY VARIABLES
793	!
794	!
795	DT4=2.*DTSTEP
796	RDT=1./DT4
797	!
798	DO I=its,ite
799	HGAMT(I)=0.
800	HGAMQ(I)=0.
801	WSCALE(I)=0.
802	KPBL(I)=KL
803	PBL(I)=ZQ(I,KL)
804	KPBL0(I)=KL
805	PBL0(I)=ZQ(I,KL)
806	PBLFLG(I)=.TRUE.
807	SFCFLG(I)=.TRUE.
808	IF(BR(I).GT.0.0)SFCFLG(I)=.FALSE.
809	ZL1(I)=ZA(I,KL)
810	THERMAL(I)=THVX(I,KL)
811	ENDDO
812
813	! COMPUTE THE FIRST GUESS OF PBL HEIGHT
814
815	DO I=its,ite
816	STABLE(I)=.FALSE.
817	BRUP(I)=BR(I)
818	ENDDO
819	DO K=KLM,KLPBL,-1
820	DO I=its,ite
821	IF(.NOT.STABLE(I))THEN
822	BRDN(I)=BRUP(I)
823	SPDK2=MAX(UX(I,K)2+VX(I,K)2,1.)
824	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
825	KPBL(I)=K
826	STABLE(I)=BRUP(I).GT.BRCR
827	ENDIF
828	ENDDO
829	ENDDO
830	!
831	DO I=its,ite
832	K=KPBL(I)
833	IF(BRDN(I).GE.BRCR)THEN
834	BRINT=0.
835	ELSEIF(BRUP(I).LE.BRCR)THEN
836	BRINT=1.
837	ELSE
838	BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
839	ENDIF
840	PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
841	IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
842	ENDDO
843	!
844	DO I=its,ite
845	FM=GZ1OZ0(I)-PSIM(I)
846	FH=GZ1OZ0(I)-PSIH(I)
847	HOL(I)=MAX(BR(I)FMFM/FH,RIMIN)
848	IF(SFCFLG(I))THEN
849	HOL(I)=MIN(HOL(I),-ZFMIN)
850	ELSE
851	HOL(I)=MAX(HOL(I),ZFMIN)
852	ENDIF
853	!
854	HOL1=HOL(I)PBL(I)/ZL1(I)SFCFRAC
855	HOL(I)=-HOL(I)*PBL(I)/ZL1(I)
856	IF(SFCFLG(I))THEN
857	PHIM(I)=(1.-APHI16HOL1)*(-1./4.)
858	PHIH(I)=(1.-APHI16HOL1)*(-1./2.)
859	ELSE
860	PHIM(I)=(1.+APHI5*HOL1)
861	PHIH(I)=PHIM(I)
862	ENDIF
863	WSCALE(I)=UST(I)/PHIM(I)
864	WSCALE(I)=MIN(WSCALE(I),UST(I)*APHI16)
865	WSCALE(I)=MAX(WSCALE(I),UST(I)/APHI5)
866	ENDDO
867
868	! COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
869	! UNDER UNSTABLE CONDITIONS
870
871	DO I=its,ite
872	IF(SFCFLG(I))THEN
873	GAMFAC=CFAC/RHOX(I)/WSCALE(I)
874	HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)
875	HGAMQ(I)=MIN(GAMFAC*QFX(I),GAMCRQ)
876	IF((XLAND(I)-1.5).GE.0)HGAMQ(I)=0.
877	VPERT=HGAMT(I)+EP1THX(I,KL)HGAMQ(I)
878	VPERT=MIN(VPERT,GAMCRT)
879	THERMAL(I)=THERMAL(I)+MAX(VPERT,0.)
880	HGAMT(I)=MAX(HGAMT(I),0.0)
881	HGAMQ(I)=MAX(HGAMQ(I),0.0)
882	ELSE
883	PBLFLG(I)=.FALSE.
884	ENDIF
885	ENDDO
886	!
887	DO I=its,ite
888	IF(PBLFLG(I))THEN
889	KPBL(I)=KL
890	PBL(I)=ZQ(I,KL)
891	ENDIF
892	ENDDO
893	!
894	! ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
895	!
896	DO I=its,ite
897	IF(PBLFLG(I))THEN
898	STABLE(I)=.FALSE.
899	BRUP(I)=BR(I)
900	ENDIF
901	ENDDO
902	DO K=KLM,KLPBL,-1
903	DO I=its,ite
904	IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
905	BRDN(I)=BRUP(I)
906	SPDK2=MAX((UX(I,K)2+VX(I,K)2),1.)
907	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
908	KPBL(I)=K
909	STABLE(I)=BRUP(I).GT.BRCR
910	ENDIF
911	ENDDO
912	ENDDO
913	!
914	DO I=its,ite
915	IF(PBLFLG(I))THEN
916	K=KPBL(I)
917	IF(BRDN(I).GE.BRCR)THEN
918	BRINT=0.
919	ELSEIF(BRUP(I).LE.BRCR)THEN
920	BRINT=1.
921	ELSE
922	BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
923	ENDIF
924	PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
925	IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
926	IF(KPBL(I).LE.1)PBLFLG(I)=.FALSE.
927	ENDIF
928	ENDDO
929	!
930	! DIAGNOSTIC PBL HEIGHT WITH BRCR EFFECTIVELY ZERO (PBL0)
931	!
932	DO I=its,ite
933	IF(PBLFLG(I))THEN
934	STABLE(I)=.FALSE.
935	BRUP(I)=BR(I)
936	ENDIF
937	ENDDO
938	DO K=KLM,KLPBL,-1
939	DO I=its,ite
940	IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
941	BRDN(I)=BRUP(I)
942	SPDK2=MAX((UX(I,K)2+VX(I,K)2),1.)
943	BRUP(I)=(THVX(I,K)-THERMAL(I))(GZA(I,K)/THVX(I,KL))/SPDK2
944	KPBL0(I)=K
945	STABLE(I)=BRUP(I).GT.0.0
946	ENDIF
947
948	ENDDO
949	ENDDO
950	!
951	DO I=its,ite
952	IF(PBLFLG(I))THEN
953	K=KPBL0(I)
954	IF(BRDN(I).GE.0.0)THEN
955	BRINT=0.
956	ELSEIF(BRUP(I).LE.0.0)THEN
957	BRINT=1.
958	ELSE
959	BRINT=(0.0-BRDN(I))/(BRUP(I)-BRDN(I))
960	ENDIF
961	PBL0(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
962	IF(PBL0(I).LT.ZQ(I,KPBL0(I)+1))KPBL0(I)=KPBL0(I)+1
963	IF(KPBL0(I).LE.1)PBLFLG(I)=.FALSE.
964	ENDIF
965	ENDDO
966
967	!
968	! COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
969	!
970	DO K=kte,KLPBL,-1
971	DO I=its,ite
972	IF(KPBL(I).LT.K)THEN
973	PRNUM=(PHIH(I)/PHIM(I)+CFACKARMANSFCFRAC)
974	PRNUM=MIN(PRNUM,PRMAX)
975	PRNUM=MAX(PRNUM,PRMIN)
976	ZFAC=MAX((1.-(ZQ(I,K)-ZL1(I))/(PBL(I)-ZL1(I))),ZFMIN)
977	XKZO=CKZ*DZA(I,K-1)
978	XKZM(I,K)=XKZO+WSCALE(I)KARMANZQ(I,K)ZFAC*PFAC
979	XKZH(I,K)=XKZM(I,K)/PRNUM
980	XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
981	XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
982	XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
983	XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
984	ENDIF
985	ENDDO
986	ENDDO
987	!
988	! COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
989	!
990	DO K=kts+1,kte
991	DO I=its,ite
992	XKZO=CKZ*DZA(I,K-1)
993	IF(K.LE.KPBL(I))THEN
994	SS=((UX(I,K-1)-UX(I,K))*(UX(I,K-1)-UX(I,K))+(VX(I,K-1)- &
995	VX(I,K))(VX(I,K-1)-VX(I,K)))/(DZA(I,K-1)DZA(I,K-1))+ &
996	1.E-9
997	RI=GOVRTH(I)(THVX(I,K-1)-THVX(I,K))/(SSDZA(I,K-1))
998	IF(IMVDIF.EQ.1)THEN
999	IF((QCX(I,K)+QIX(I,K)).GT.0.01E-3.AND.(QCX(I,K-1)+ &
1000	QIX(I,K-1)).GT.0.01E-3)THEN
1001	! IN CLOUD
1002	QMEAN=0.5*(QX(I,K)+QX(I,K-1))
1003	TMEAN=0.5*(SCR3(I,K)+SCR3(I,K-1))
1004	ALPH=XLV*QMEAN/R/TMEAN
1005	CHI=XLVXLVQMEAN/CP/RV/TMEAN/TMEAN
1006	RI=(1.+ALPH)(RI-GG/SS/TMEAN/CP*((CHI-ALPH)/(1.+CHI)))
1007	ENDIF
1008	ENDIF
1009	ZK=KARMAN*ZQ(I,K)
1010	RL2=(ZKRLAM/(RLAM+ZK))*2
1011	DK=RL2*SQRT(SS)
1012	IF(RI.LT.0.)THEN
1013	! UNSTABLE REGIME
1014	SRI=SQRT(-RI)
1015	XKZM(I,K)=XKZO+DK(1+8.(-RI)/(1+1.746*SRI))
1016	XKZH(I,K)=XKZO+DK(1+8.(-RI)/(1+1.286*SRI))
1017	ELSE
1018	! STABLE REGIME
1019	XKZH(I,K)=XKZO+DK/(1+5.RI)*2
1020	PRNUM=1.0+2.1*RI
1021	PRNUM=MIN(PRNUM,PRMAX)
1022	XKZM(I,K)=(XKZH(I,K)-XKZO)*PRNUM+XKZO
1023	ENDIF
1024	!
1025	XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
1026	XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
1027	XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
1028	XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
1029	ENDIF
1030	!
1031	ENDDO
1032	ENDDO
1033
1034	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
1035
1036	DO I=its,ite
1037	DO K=kts,kte
1038	AU(i,k)=0.
1039	AL(i,k)=0.
1040	AD(i,k)=0.
1041	A1(i,k)=0.
1042	A2(i,k)=0.
1043	ENDDO
1044	ENDDO
1045
1046	DO I=its,ite
1047	AD(I,1)=1.
1048	A1(I,1)=SCR3(I,KL)+HFX(I)/(RHOX(I)CPM(I))/ZQ(I,KL)DT4
1049	A2(I,1)=QX(I,KL)+QFX(I)/(RHOX(I))/ZQ(I,KL)*DT4
1050	ENDDO
1051	!
1052	DO K=kte,kts+1,-1
1053	KK=kme-K
1054	DO I=its,ite
1055	DTODSD=DT4/dz8w2d(I,K)
1056	DTODSU=DT4/dz8w2d(I,K-1)
1057	DSIG=z2d(I,K)-z2d(I,K-1)
1058	DSIG=-DSIG
1059	RDZ=1./DZA(I,K-1)
1060	IF(PBLFLG(I).AND.KPBL(I).LT.K)THEN
1061	DSDZT=DSIGXKZH(I,K)RDZ*(G/CP-HGAMT(I)/PBL(I))
1062	DSDZQ=DSIGXKZH(I,K)RDZ*(-HGAMQ(I)/PBL(I))
1063	A2(I,KK)=A2(I,KK)+DTODSD*DSDZQ
1064	A2(I,KK+1)=QX(I,K-1)-DTODSU*DSDZQ
1065	ELSE
1066	DSDZT=DSIGXKZH(I,K)RDZ*(G/CP)
1067	A2(I,KK+1)=QX(I,K-1)
1068	ENDIF
1069	DSDZ2=DSIGXKZH(I,K)RDZ*RDZ
1070	AU(I,KK)=-DTODSD*DSDZ2
1071	AL(I,KK)=-DTODSU*DSDZ2
1072	AD(I,KK)=AD(I,KK)-AU(I,KK)
1073	AD(I,KK+1)=1.-AL(I,KK)
1074	A1(I,KK)=A1(I,KK)+DTODSD*DSDZT
1075	A1(I,KK+1)=SCR3(I,K-1)-DTODSU*DSDZT
1076	ENDDO
1077	ENDDO
1078
1079	! SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
1080
1081	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1082	its,ite,kts,kte )
1083
1084	! RECOVER TENDENCIES OF HEAT AND MOISTURE
1085
1086	DO K=kte,kts,-1
1087	KK=kme-K
1088	DO I=its,ite
1089	TTEND=(A1(I,KK)-SCR3(I,K))*RDT
1090	QTEND=(A2(I,KK)-QX(I,K))*RDT
1091	TTNP(I,K)=TTNP(I,K)+TTEND
1092	QTNP(I,K)=QTNP(I,K)+QTEND
1093	ENDDO
1094	ENDDO
1095
1096	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
1097
1098	DO I=its,ite
1099	DO K=kts,kte
1100	AU(i,k)=0.
1101	AL(i,k)=0.
1102	AD(i,k)=0.
1103	A1(i,k)=0.
1104	A2(i,k)=0.
1105	ENDDO
1106	ENDDO
1107
1108	DO I=its,ite
1109	AD(I,1)=1.
1110	A1(I,1)=UX(I,KL)-UX(I,KL)/WSPD1(I)UST(I)UST(I)/ZQ(I,KL) &
1111	DT4(WSPD1(I)/WSPD(I))**2
1112	A2(I,1)=VX(I,KL)-VX(I,KL)/WSPD1(I)UST(I)UST(I)/ZQ(I,KL) &
1113	DT4(WSPD1(I)/WSPD(I))**2
1114	ENDDO
1115	!
1116	DO K=kte,kts+1,-1
1117	KK=kme-K
1118	DO I=its,ite
1119	DTODSD=DT4/dz8w2d(I,K)
1120	DTODSU=DT4/dz8w2d(I,K-1)
1121	DSIG=z2d(I,K)-z2d(I,K-1)
1122	DSIG=-DSIG
1123	RDZ=1./DZA(I,K-1)
1124	DSDZ2=DSIGXKZM(I,K)RDZ*RDZ
1125	AU(I,KK)=-DTODSD*DSDZ2
1126	AL(I,KK)=-DTODSU*DSDZ2
1127	AD(I,KK)=AD(I,KK)-AU(I,KK)
1128	AD(I,KK+1)=1.-AL(I,KK)
1129	A1(I,KK+1)=UX(I,K-1)
1130	A2(I,KK+1)=VX(I,K-1)
1131	ENDDO
1132	ENDDO
1133
1134	! SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
1135
1136	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1137	its,ite,kts,kte )
1138
1139	! RECOVER TENDENCIES OF MOMENTUM
1140
1141	DO K=kte,kts,-1
1142	KK=kme-K
1143	DO I=its,ite
1144	UTEND=(A1(I,KK)-UX(I,K))*RDT
1145	VTEND=(A2(I,KK)-VX(I,K))*RDT
1146	UTNP(I,K)=UTNP(I,K)+UTEND
1147	VTNP(I,K)=VTNP(I,K)+VTEND
1148	ENDDO
1149	ENDDO
1150
1151	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR CLOUD
1152
1153	DO I=its,ite
1154	DO K=kts,kte
1155	AU(i,k)=0.
1156	AL(i,k)=0.
1157	AD(i,k)=0.
1158	A1(i,k)=0.
1159	A2(i,k)=0.
1160	ENDDO
1161	ENDDO
1162
1163	! IF(IMOIST.EQ.1)GOTO 690
1164	DO I=its,ite
1165	AD(I,1)=1.
1166	A1(I,1)=QCX(I,KL)
1167	A2(I,1)=QIX(I,KL)
1168	ENDDO
1169	!
1170	DO K=kte,kts+1,-1
1171	KK=kme-K
1172	DO I=its,ite
1173	DTODSD=DT4/dz8w2d(I,K)
1174	DTODSU=DT4/dz8w2d(I,K-1)
1175	DSIG=z2d(I,K)-z2d(I,K-1)
1176	DSIG=-DSIG
1177	RDZ=1./DZA(I,K-1)
1178	A1(I,KK+1)=QCX(I,K-1)
1179	A2(I,KK+1)=QIX(I,K-1)
1180	DSDZ2=DSIGXKZH(I,K)RDZ*RDZ
1181	AU(I,KK)=-DTODSD*DSDZ2
1182	AL(I,KK)=-DTODSU*DSDZ2
1183	AD(I,KK)=AD(I,KK)-AU(I,KK)
1184	AD(I,KK+1)=1.-AL(I,KK)
1185	ENDDO
1186	ENDDO
1187
1188	! SOLVE TRIDIAGONAL PROBLEM FOR CLOUD
1189
1190	CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1191	its,ite,kts,kte )
1192	!
1193	DO K=kte,kts,-1
1194	KK=kme-K
1195	DO I=its,ite
1196	QCTEND=(A1(I,KK)-QCX(I,K))*RDT
1197	QITEND=(A2(I,KK)-QIX(I,K))*RDT
1198	QCTNP(I,K)=QCTNP(I,K)+QCTEND
1199	QITNP(I,K)=QITNP(I,K)+QITEND
1200	ENDDO
1201	ENDDO
1202	!
1203	!---- END OF VERTICAL DIFFUSION
1204	!
1205	690 CONTINUE
1206	!
1207	!-----CALCULATION OF NEW VALUES DUE TO VERTICAL EXCHANGE PROCESSES IS
1208	! COMPLETED. THE FINAL STEP IS TO ADD THE TENDENCIES CALCULATED
1209	! IN HIRPBL TO THOSE OF MM4.
1210
1211	DO 820 K=kts,kte
1212	NK=kme-K
1213	DO 820 I=its,ite
1214	U2DTEN(I,NK)=UTNP(I,K)
1215	V2DTEN(I,NK)=VTNP(I,K)
1216	820 CONTINUE
1217	!
1218	! IF(J.EQ.1.AND.IN.GT.1)GOTO 860
1219	!SUE JBGN=3
1220	!SUE JEND=JLXM-1
1221
1222	! change when nest
1223	!SUE JBGN=2
1224	!SUE JEND=JLXM
1225
1226	JBGN=jts
1227	JEND=jte
1228	IBGN=its
1229	IEND=ite
1230
1231	! IF(J.LT.JBGN.OR.J.GT.JEND)GOTO 860
1232	!SUE IBGN=3
1233	!SUE IEND=ILXM-1
1234
1235	! change when nest
1236	!SUE IBGN=2
1237	!SUE IEND=ILXM
1238
1239	DO 830 K=kts,kte
1240	NK=kme-K
1241	DO 830 I=IBGN,IEND
1242	T2DTEN(I,NK)=TTNP(I,K)
1243	830 CONTINUE
1244	!
1245	! IF(IDRY.EQ.1)GOTO 860
1246	DO 840 K=kts,kte
1247	NK=kme-K
1248	DO 840 I=IBGN,IEND
1249	QV2DTEN(I,NK)=QTNP(I,K)
1250	840 CONTINUE
1251
1252	! IF(IMOIST.EQ.1)GOTO 860
1253	DO 850 K=kts,kte
1254	NK=kme-K
1255	DO 850 I=IBGN,IEND
1256	QC2DTEN(I,NK)=QCTNP(I,K)
1257	850 CONTINUE
1258
1259	IF(flag_QI .AND. PRESENT( QI2DTEN ) ) THEN
1260	DO K=kts,kte
1261	NK=kme-K
1262	DO I=IBGN,IEND
1263	QI2DTEN(I,NK)=QITNP(I,K)
1264	ENDDO
1265	ENDDO
1266	ENDIF
1267
1268	860 CONTINUE
1269	!
1270	!-----APPLY ASSELIN FILTER TO TGD FOR LARGE TIME STEP:
1271	!
1272	! DO 885 I=its,ite
1273	! TSK(I)=TSK(I)(PS(I)/100.)*ROVCP
1274	! 885 CONTINUE
1275	!
1276	940 CONTINUE
1277	!
1278	! KPBL IS NEEDED FOR THE FDDA, AND SINCE THERE IS NO LONGER JUST ONE
1279	! LARGE "J LOOP" IT MUST BE STORED AS (I,J)...
1280	!
1281	! USE NEW DIAGNOSED PBL DEPTH (CALCULATED WITH brcr=0.0)
1282	! PBL IS USED FOR OUTPUT AND NEXT-TIME-STEP BELJAARS CALC IN SFCLAY
1283	DO 950 I=its,ite
1284	KPBL1D(I)=KPBL0(I)
1285	PBL(I)=PBL0(I)
1286	950 CONTINUE
1287
1288	END SUBROUTINE MRF2D
1289
1290	!================================================================
1291	SUBROUTINE TRIDI2(CL,CM,CU,R1,R2,AU,A1,A2, &
1292	its,ite,kts,kte )
1293	!----------------------------------------------------------------
1294	IMPLICIT NONE
1295	!----------------------------------------------------------------
1296
1297	INTEGER, INTENT(IN ) :: its,ite, kts,kte
1298
1299	REAL, DIMENSION( its:ite, kts+1:kte+1 ) , &
1300	INTENT(IN ) :: CL
1301
1302	REAL, DIMENSION( its:ite, kts:kte ) , &
1303	INTENT(IN ) :: CM, &
1304	R1, &
1305	R2
1306	REAL, DIMENSION( its:ite, kts:kte ) , &
1307	INTENT(INOUT) :: AU, &
1308	CU, &
1309	A1, &
1310	A2
1311
1312	REAL :: FK
1313	INTEGER :: I,K,L,N
1314
1315	!----------------------------------------------------------------
1316
1317	L=ite
1318	N=kte
1319
1320	DO I=its,L
1321	FK=1./CM(I,1)
1322	AU(I,1)=FK*CU(I,1)
1323	A1(I,1)=FK*R1(I,1)
1324	A2(I,1)=FK*R2(I,1)
1325	ENDDO
1326	DO K=2,N-1
1327	DO I=its,L
1328	FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1329	AU(I,K)=FK*CU(I,K)
1330	A1(I,K)=FK(R1(I,K)-CL(I,K)A1(I,K-1))
1331	A2(I,K)=FK(R2(I,K)-CL(I,K)A2(I,K-1))
1332	ENDDO
1333	ENDDO
1334	DO I=its,L
1335	FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1336	A1(I,N)=FK(R1(I,N)-CL(I,N)A1(I,N-1))
1337	A2(I,N)=FK(R2(I,N)-CL(I,N)A2(I,N-1))
1338
1339	ENDDO
1340	DO K=N-1,kts,-1
1341	DO I=its,L
1342	A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
1343	A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
1344	ENDDO
1345	ENDDO
1346
1347	END SUBROUTINE TRIDI2
1348
1349	!===================================================================
1350	SUBROUTINE mrfinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1351	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
1352	restart, allowed_to_read , &
1353	ids, ide, jds, jde, kds, kde, &
1354	ims, ime, jms, jme, kms, kme, &
1355	its, ite, jts, jte, kts, kte )
1356	!-------------------------------------------------------------------
1357	IMPLICIT NONE
1358	!-------------------------------------------------------------------
1359	LOGICAL , INTENT(IN) :: restart , allowed_to_read
1360	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1361	ims, ime, jms, jme, kms, kme, &
1362	its, ite, jts, jte, kts, kte
1363	INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
1364
1365	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
1366	RUBLTEN, &
1367	RVBLTEN, &
1368	RTHBLTEN, &
1369	RQVBLTEN, &
1370	RQCBLTEN, &
1371	RQIBLTEN
1372	INTEGER :: i, j, k, itf, jtf, ktf
1373
1374	jtf=min0(jte,jde-1)
1375	ktf=min0(kte,kde-1)
1376	itf=min0(ite,ide-1)
1377
1378	IF(.not.restart)THEN
1379	DO j=jts,jtf
1380	DO k=kts,ktf
1381	DO i=its,itf
1382	RUBLTEN(i,k,j)=0.
1383	RVBLTEN(i,k,j)=0.
1384	RTHBLTEN(i,k,j)=0.
1385	RQVBLTEN(i,k,j)=0.
1386	RQCBLTEN(i,k,j)=0.
1387	ENDDO
1388	ENDDO
1389	ENDDO
1390	ENDIF
1391
1392	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
1393	DO j=jts,jtf
1394	DO k=kts,ktf
1395	DO i=its,itf
1396	RQIBLTEN(i,k,j)=0.
1397	ENDDO
1398	ENDDO
1399	ENDDO
1400	ENDIF
1401
1402	END SUBROUTINE mrfinit
1403
1404	!-------------------------------------------------------------------
1405
1406	END MODULE module_bl_mrf
1407

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

Download in other formats:

Original Format