Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_mrf.F @ 3567

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

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

Download in other formats:

Original Format