Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_myjpbl.F @ 3431

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

Line
1	!-----------------------------------------------------------------------
2	!
3	MODULE MODULE_BL_MYJPBL
4	!
5	!-----------------------------------------------------------------------
6	!
7	USE MODULE_MODEL_CONSTANTS
8	!
9	!-----------------------------------------------------------------------
10	!
11	! REFERENCES: Janjic (2002), NCEP Office Note 437
12	! Mellor and Yamada (1982), Rev. Geophys. Space Phys.
13	!
14	! ABSTRACT:
15	! MYJ UPDATES THE TURBULENT KINETIC ENERGY WITH THE PRODUCTION/
16	! DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
17	! (USING AN IMPLICIT FORMULATION) FROM MELLOR-YAMADA
18	! LEVEL 2.5 AS EXTENDED BY JANJIC. EXCHANGE COEFFICIENTS FOR
19	! THE SURFACE AND FOR ALL LAYER INTERFACES ARE COMPUTED FROM
20	! MONIN-OBUKHOV THEORY.
21	! THE TURBULENT VERTICAL EXCHANGE IS THEN EXECUTED.
22	!
23	!-----------------------------------------------------------------------
24	!
25	INTEGER :: ITRMX=5 ! Iteration count for mixing length computation
26	!
27	! REAL,PARAMETER :: G=9.81,PI=3.1415926,R_D=287.04,R_V=461.6 &
28	! & ,VKARMAN=0.4
29	REAL,PARAMETER :: PI=3.1415926,VKARMAN=0.4
30	! REAL,PARAMETER :: CP=7.*R_D/2.
31	REAL,PARAMETER :: CAPA=R_D/CP
32	REAL,PARAMETER :: RLIVWV=XLS/XLV,ELOCP=2.72E6/CP
33	REAL,PARAMETER :: EPS1=1.E-12,EPS2=0.
34	REAL,PARAMETER :: EPSL=0.32,EPSRU=1.E-7,EPSRS=1.E-7 &
35	& ,EPSTRB=1.E-24
36	REAL,PARAMETER :: EPSA=1.E-8,EPSIT=1.E-4,EPSU2=1.E-4,EPSUST=0.07 &
37	& ,FH=1.01
38	REAL,PARAMETER :: ALPH=0.30,BETA=1./273.,EL0MAX=1000.,EL0MIN=1. &
39	& ,ELFC=0.23*0.5,GAM1=0.2222222222222222222 &
40	& ,PRT=1.
41	REAL,PARAMETER :: A1=0.659888514560862645 &
42	& ,A2x=0.6574209922667784586 &
43	& ,B1=11.87799326209552761 &
44	& ,B2=7.226971804046074028 &
45	& ,C1=0.000830955950095854396
46	REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
47	REAL,PARAMETER :: ELZ0=0.,ESQ=5.0,EXCM=0.001 &
48	& ,FHNEU=0.8,GLKBR=10.,GLKBS=30. &
49	& ,QVISC=2.1E-5,RFC=0.191,RIC=0.505,SMALL=0.35 &
50	& ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5 &
51	& ,USTC=0.7,USTR=0.225,VISC=1.5E-5 &
52	& ,WOLD=0.15,WWST=1.2,ZTMAX=1.,ZTFC=1.,ZTMIN=-5.
53	!
54	REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
55	!
56	REAL,PARAMETER :: BTG=BETAG,CZIV=SMALLGLKBS &
57	! & ,EP_1=R_V/R_D-1.,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
58	& ,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
59	& ,RB1=1./B1,RTVISC=1./TVISC,RVISC=1./VISC &
60	& ,ZQRZT=SQSC/SQPR
61	!
62	REAL,PARAMETER :: ADNH= 9.A1A2xA2x(12.A1+3.B2)BTGBTG &
63	& ,ADNM=18.A1A1A2x(B2-3.A2x)BTG &
64	& ,ANMH=-9.A1A2xA2xBTG*BTG &
65	& ,ANMM=-3.A1A2x(3.A2x+3.B2C1+18.A1C1-B2) &
66	& *BTG &
67	& ,BDNH= 3.A2x(7.A1+B2)BTG &
68	& ,BDNM= 6.A1A1 &
69	& ,BEQH= A2xB1BTG+3.A2x(7.A1+B2)BTG &
70	& ,BEQM=-A1B1(1.-3.C1)+6.A1*A1 &
71	& ,BNMH=-A2x*BTG &
72	& ,BNMM=A1(1.-3.C1) &
73	& ,BSHH=9.A1A2xA2xBTG &
74	& ,BSHM=18.A1A1A2xC1 &
75	& ,BSMH=-3.A1A2x(3.A2x+3.B2C1+12.A1C1-B2) &
76	& *BTG &
77	& ,CESH=A2x &
78	& ,CESM=A1(1.-3.C1) &
79	& ,CNV=EP_1*G/BTG &
80	& ,ELFCS=VKARMAN*BTG &
81	& ,FZQ1=RTVISCQVISCZQRZT &
82	& ,FZQ2=RTVISCQVISCZQRZT &
83	& ,FZT1=RVISC TVISCSQPR &
84	& ,FZT2=CZIVGRRSTVISC*SQPR &
85	& ,FZU1=CZIV*VISC &
86	& ,PIHF=0.5*PI &
87	& ,RFAC=RIC/(FHNEURFCRFC) &
88	& ,RQVISC=1./QVISC &
89	& ,RRIC=1./RIC &
90	& ,USTFC=0.018/G &
91	& ,WNEW=1.-WOLD &
92	& ,WWST2=WWST*WWST
93	!
94	!-----------------------------------------------------------------------
95	!* FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)2
96	!-----------------------------------------------------------------------
97	!
98	REAL,PARAMETER :: AEQH=9.A1A2xA2xB1BTGBTG &
99	& +9.A1A2xA2x(12.A1+3.B2)BTGBTG &
100	& ,AEQM=3.A1A2xB1(3.A2x+3.B2C1+18.A1*C1-B2)&
101	& BTG+18.A1A1A2x(B2-3.A2x)*BTG
102	!
103	!-----------------------------------------------------------------------
104	!*** FORBIDDEN TURBULENCE AREA
105	!-----------------------------------------------------------------------
106	!
107	REAL,PARAMETER :: REQU=-AEQH/AEQM &
108	& ,EPSGH=1.E-9,EPSGM=REQU*EPSGH
109	!
110	!-----------------------------------------------------------------------
111	!*** NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
112	!-----------------------------------------------------------------------
113	!
114	REAL,PARAMETER :: UBRYL=(18.REQUA1A1A2xB2C1*BTG &
115	& +9.A1A2xA2xB2BTGBTG) &
116	& /(REQU*ADNM+ADNH) &
117	& ,UBRY=(1.+EPSRS)UBRYL,UBRY3=3.UBRY
118	!
119	REAL,PARAMETER :: AUBH=27.A1A2xA2xB2BTGBTG-ADNH*UBRY3 &
120	& ,AUBM=54.A1A1A2xB2C1BTG -ADNM*UBRY3 &
121	& ,BUBH=(9.A1A2x+3.A2xB2)BTG-BDNHUBRY3 &
122	& ,BUBM=18.A1A1C1 -BDNMUBRY3 &
123	& ,CUBR=1. - UBRY3 &
124	& ,RCUBR=1./CUBR
125	!
126	!-----------------------------------------------------------------------
127	!
128	CONTAINS
129	!
130	!----------------------------------------------------------------------
131	SUBROUTINE MYJPBL(DT,STEPBL,HT,DZ &
132	& ,PMID,PINT,TH,T,EXNER,QV,CWM,U,V,RHO &
133	& ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 &
134	& ,LOWLYR,XLAND,SICE,SNOW &
135	& ,TKE_MYJ,EXCH_H,USTAR,ZNT,EL_MYJ,PBLH,KPBL,CT &
136	& ,AKHS,AKMS,ELFLX &
137	& ,RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN &
138	& ,IDS,IDE,JDS,JDE,KDS,KDE &
139	& ,IMS,IME,JMS,JME,KMS,KME &
140	& ,ITS,ITE,JTS,JTE,KTS,KTE)
141	!----------------------------------------------------------------------
142	!
143	IMPLICIT NONE
144	!
145	!----------------------------------------------------------------------
146	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
147	& ,IMS,IME,JMS,JME,KMS,KME &
148	& ,ITS,ITE,JTS,JTE,KTS,KTE
149	!
150	INTEGER,INTENT(IN) :: STEPBL
151
152	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
153	!
154	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: KPBL
155	!
156	REAL,INTENT(IN) :: DT
157	!
158	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,SICE,SNOW &
159	& ,TSK,XLAND
160	!
161	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DZ &
162	& ,EXNER &
163	& ,PMID,PINT &
164	& ,QV,RHO &
165	& ,T,TH,U,V
166	!
167	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PBLH
168	!
169	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS
170	!
171	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
172	& ,INTENT(OUT) :: EL_MYJ &
173	& ,RQCBLTEN,RQVBLTEN &
174	& ,RTHBLTEN &
175	& ,RUBLTEN,RVBLTEN
176	!
177	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CT,QSFC,QZ0 &
178	& ,THZ0,USTAR &
179	& ,UZ0,VZ0,ZNT
180	!
181	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
182	& ,INTENT(INOUT) :: EXCH_H,TKE_MYJ
183	!
184	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHKLOWQ,ELFLX
185	!
186	!----------------------------------------------------------------------
187	!***
188	!*** LOCAL VARIABLES
189	!***
190	INTEGER :: I,J,K,KFLIP,LLOW,LMH,LMXL
191	!
192	INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: LPBL
193	!
194	REAL :: AKHS_DENS,AKMS_DENS,APEX,DCDT,DELTAZ,DQDT,DTDIF,DTDT &
195	& ,DTTURBL,DUDT,DVDT,EXNSFC,PSFC,PTOP,QFC1,QLOW,QOLD &
196	& ,RATIOMX,RDTTURBL,RG,RWMSK,SEAMASK,THNEW,THOLD,TX &
197	& ,ULOW,VLOW,WMSK
198	!
199	REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,TK,UK,VK
200	!
201	REAL,DIMENSION(KTS:KTE-1) :: AKHK,AKMK,EL,GH,GM
202	!
203	REAL,DIMENSION(KTS:KTE+1) :: ZHK
204	!
205	REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
206	!
207	REAL,DIMENSION(KTS:KTE,ITS:ITE) :: RHOK
208	!
209	REAL,DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: APE,THE
210	!
211	REAL,DIMENSION(ITS:ITE,KTS:KTE-1,JTS:JTE) :: AKH,AKM
212	!
213	REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
214	!
215	!*** Begin debugging
216	REAL :: ZSL_DIAG
217	INTEGER :: IMD,JMD,PRINT_DIAG
218	!*** End debugging
219	!
220	!----------------------------------------------------------------------
221	!**********************************************************************
222	!----------------------------------------------------------------------
223	!
224	!*** Begin debugging
225	IMD=(IMS+IME)/2
226	JMD=(JMS+JME)/2
227	!*** End debugging
228	!
229	!*** MAKE PREPARATIONS
230	!
231	!----------------------------------------------------------------------
232	DTTURBL=DT*STEPBL
233	RDTTURBL=1./DTTURBL
234	DTDIF=DTTURBL
235	RG=1./G
236	!
237	DO J=JTS,JTE
238	DO K=KTS,KTE-1
239	DO I=ITS,ITE
240	AKM(I,K,J)=0.
241	ENDDO
242	ENDDO
243	ENDDO
244	!
245	DO J=JTS,JTE
246	DO K=KTS,KTE+1
247	DO I=ITS,ITE
248	ZINT(I,K,J)=0.
249	ENDDO
250	ENDDO
251	ENDDO
252	!
253	DO J=JTS,JTE
254	DO I=ITS,ITE
255	ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer
256	!
257	!!!!!!!!!
258	!!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
259	!!!!!!!!!
260	!!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer
261	!!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer
262	!
263	ENDDO
264	ENDDO
265	!
266	DO J=JTS,JTE
267	DO K=KTE,KTS,-1
268	KFLIP=KTE+1-K
269	DO I=ITS,ITE
270	ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
271	APEX=1./EXNER(I,K,J)
272	APE(I,K,J)=APEX
273	TX=T(I,K,J)
274	THE(I,K,J)=(CWM(I,K,J)(-ELOCP/TX)+1.)TH(I,K,J)
275	ENDDO
276	ENDDO
277	ENDDO
278	!
279	EL_MYJ = 0.
280	!
281	!----------------------------------------------------------------------
282	setup_integration: DO J=JTS,JTE
283	!----------------------------------------------------------------------
284	!
285	DO I=ITS,ITE
286	!
287	!*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
288	!
289	LMH=KTE-LOWLYR(I,J)+1
290	!
291	PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
292	PSFC=PINT(I,LOWLYR(I,J),J)
293	!
294	!*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
295	!
296	SEAMASK=XLAND(I,J)-1.
297	!
298	!*** FILL 1-D VERTICAL ARRAYS
299	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
300	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
301	!
302	DO K=KTE,KTS,-1
303	KFLIP=KTE+1-K
304	TK(K)=T(I,KFLIP,J)
305	THEK(K)=THE(I,KFLIP,J)
306	RATIOMX=QV(I,KFLIP,J)
307	QK(K)=RATIOMX/(1.+RATIOMX)
308	CWMK(K)=CWM(I,KFLIP,J)
309	PK(K)=PMID(I,KFLIP,J)
310	UK(K)=U(I,KFLIP,J)
311	VK(K)=V(I,KFLIP,J)
312	!
313	!*** TKE=0.5(q2) ==> q2=2.TKE
314	!
315	Q2K(K)=2.*TKE_MYJ(I,KFLIP,J)
316	!
317	!*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
318	!
319	ZHK(K)=ZINT(I,K,J)
320	!
321	ENDDO
322	ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer
323	!
324	!*** Begin debugging
325	! IF(I==IMD.AND.J==JMD)THEN
326	! PRINT_DIAG=1
327	! ELSE
328	! PRINT_DIAG=0
329	! ENDIF
330	! IF(I==227.AND.J==363)PRINT_DIAG=2
331	!*** End debugging
332	!
333	!----------------------------------------------------------------------
334	!***
335	!*** FIND THE MIXING LENGTH
336	!***
337	CALL MIXLEN(LMH,UK,VK,TK,THEK,QK,CWMK &
338	& ,Q2K,ZHK,GM,GH,EL &
339	& ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J) &
340	& ,IDS,IDE,JDS,JDE,KDS,KDE &
341	& ,IMS,IME,JMS,JME,KMS,KME &
342	& ,ITS,ITE,JTS,JTE,KTS,KTE)
343	!
344	!----------------------------------------------------------------------
345	!***
346	!*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
347	!*** THE TURBULENT KINETIC ENERGY
348	!***
349	!
350	CALL PRODQ2(LMH,DTTURBL,USTAR(I,J),GM,GH,EL,Q2K &
351	& ,IDS,IDE,JDS,JDE,KDS,KDE &
352	& ,IMS,IME,JMS,JME,KMS,KME &
353	& ,ITS,ITE,JTS,JTE,KTS,KTE)
354	!
355	!----------------------------------------------------------------------
356	!*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
357	!----------------------------------------------------------------------
358	!
359	KPBL(I,J)=KTE-LPBL(I,J)+1
360	!
361	!----------------------------------------------------------------------
362	!***
363	!*** FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
364	!***
365	CALL DIFCOF(LMH,LMXL,GM,GH,EL,TK,Q2K,ZHK,AKMK,AKHK &
366	& ,IDS,IDE,JDS,JDE,KDS,KDE &
367	& ,IMS,IME,JMS,JME,KMS,KME &
368	& ,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
369	!
370	!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
371	!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. COUNTING
372	!*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
373	!*** ARE DEFINED ON THE TOPS OF THE LAYERS KTS TO KTE-1.
374	!
375	DO K=KTS,KTE-1
376	KFLIP=KTE-K
377	AKH(I,K,J)=AKHK(K)
378	AKM(I,K,J)=AKMK(K)
379	DELTAZ=0.5*(ZHK(KFLIP)-ZHK(KFLIP+2))
380	EXCH_H(I,K,J)=AKHK(KFLIP)*DELTAZ
381	ENDDO
382	!
383	!----------------------------------------------------------------------
384	!***
385	!*** CARRY OUT THE VERTICAL DIFFUSION OF
386	!*** TURBULENT KINETIC ENERGY
387	!***
388	!
389	CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK &
390	& ,IDS,IDE,JDS,JDE,KDS,KDE &
391	& ,IMS,IME,JMS,JME,KMS,KME &
392	& ,ITS,ITE,JTS,JTE,KTS,KTE)
393	!
394	!*** SAVE THE NEW TKE AND MIXING LENGTH.
395	!
396	DO K=KTS,KTE
397	KFLIP=KTE+1-K
398	Q2K(KFLIP)=AMAX1(Q2K(KFLIP),EPSQ2)
399	TKE_MYJ(I,K,J)=0.5*Q2K(KFLIP)
400	IF(KFLIP/=KTE)EL_MYJ(I,K,J)=EL(KFLIP) ! EL IS NOT DEFINED AT KTE (ground surface)
401	ENDDO
402	!
403	ENDDO
404	!
405	!----------------------------------------------------------------------
406	ENDDO setup_integration
407	!----------------------------------------------------------------------
408	!
409	!*** CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
410	!
411	DO J=JTS,JTE
412	DO I=ITS,ITE
413	PSFC=PINT(I,LOWLYR(I,J),J)
414	THSK(I,J)=TSK(I,J)(1.E5/PSFC)*CAPA
415	ENDDO
416	ENDDO
417	!
418	!----------------------------------------------------------------------
419	!
420	!----------------------------------------------------------------------
421	main_integration: DO J=JTS,JTE
422	!----------------------------------------------------------------------
423	!
424	DO I=ITS,ITE
425	!
426	!*** FILL 1-D VERTICAL ARRAYS
427	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
428	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
429	!
430	DO K=KTE,KTS,-1
431	KFLIP=KTE+1-K
432	THEK(K)=THE(I,KFLIP,J)
433	RATIOMX=QV(I,KFLIP,J)
434	QK(K)=RATIOMX/(1.+RATIOMX)
435	CWMK(K)=CWM(I,KFLIP,J)
436	ZHK(K)=ZINT(I,K,J)
437	RHOK(K,I)=PMID(I,KFLIP,J)/(R_DT(I,KFLIP,J) &
438	& (1.+P608*QK(K)-CWMK(K)))
439	ENDDO
440	!
441	!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
442	!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. THESE COEFFICIENTS
443	!*** ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
444	!
445	DO K=KTS,KTE-1
446	AKHK(K)=AKH(I,K,J)0.5(RHOK(K,I)+RHOK(K+1,I))
447	ENDDO
448	!
449	ZHK(KTE+1)=ZINT(I,KTE+1,J)
450	!
451	SEAMASK=XLAND(I,J)-1.
452	THZ0(I,J)=(1.-SEAMASK)THSK(I,J)+SEAMASKTHZ0(I,J)
453	!
454	LLOW=LOWLYR(I,J)
455	AKHS_DENS=AKHS(I,J)*RHOK(KTE+1-LLOW,I)
456	!
457	IF(SEAMASK<0.5)THEN
458	QFC1=XLVCHKLOWQ(I,J)AKHS_DENS
459	!
460	IF(SNOW(I,J)>0..OR.SICE(I,J)>0.5)THEN
461	QFC1=QFC1*RLIVWV
462	ENDIF
463	!
464	IF(QFC1>0.)THEN
465	QLOW=QK(KTE+1-LLOW)
466	QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
467	ENDIF
468	!
469	ELSE
470	PSFC=PINT(I,LOWLYR(I,J),J)
471	EXNSFC=(1.E5/PSFC)**CAPA
472	QSFC(I,J)=PQ0SEA/PSFC &
473	& EXP(A2(THSK(I,J)-A3EXNSFC)/(THSK(I,J)-A4EXNSFC))
474	ENDIF
475	!
476	QZ0 (I,J)=(1.-SEAMASK)QSFC(I,J)+SEAMASKQZ0 (I,J)
477	!
478	!*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
479	!
480	LMH=KTE-LOWLYR(I,J)+1
481	!
482	!----------------------------------------------------------------------
483	!*** CARRY OUT THE VERTICAL DIFFUSION OF
484	!*** TEMPERATURE AND WATER VAPOR
485	!----------------------------------------------------------------------
486	!
487	CALL VDIFH(DTDIF,LMH,THZ0(I,J),QZ0(I,J) &
488	& ,AKHS_DENS,CHKLOWQ(I,J),CT(I,J) &
489	& ,THEK,QK,CWMK,AKHK,ZHK,RHOK(KTS,I) &
490	& ,IDS,IDE,JDS,JDE,KDS,KDE &
491	& ,IMS,IME,JMS,JME,KMS,KME &
492	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
493	!----------------------------------------------------------------------
494	!***
495	!*** COMPUTE PRIMARY VARIABLE TENDENCIES
496	!***
497	DO K=KTS,KTE
498	KFLIP=KTE+1-K
499	THOLD=TH(I,K,J)
500	THNEW=THEK(KFLIP)+CWMK(KFLIP)ELOCPAPE(I,K,J)
501	DTDT=(THNEW-THOLD)*RDTTURBL
502	QOLD=QV(I,K,J)/(1.+QV(I,K,J))
503	DQDT=(QK(KFLIP)-QOLD)*RDTTURBL
504	DCDT=(CWMK(KFLIP)-CWM(I,K,J))*RDTTURBL
505	!
506	RTHBLTEN(I,K,J)=DTDT
507	RQVBLTEN(I,K,J)=DQDT/(1.-QK(KFLIP))**2
508	RQCBLTEN(I,K,J)=DCDT
509	ENDDO
510	!
511	!*** Begin debugging
512	! IF(I==IMD.AND.J==JMD)THEN
513	! PRINT_DIAG=0
514	! ELSE
515	! PRINT_DIAG=0
516	! ENDIF
517	! IF(I==227.AND.J==363)PRINT_DIAG=0
518	!*** End debugging
519	!
520	PSFC=.01*PINT(I,LOWLYR(I,J),J)
521	ZSL_DIAG=0.5*DZ(I,1,J)
522	!
523	!*** Begin debugging
524	! IF(PRINT_DIAG==1)THEN
525	!
526	! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
527	! '{turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
528	! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
529	! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
530	! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
531	! '{turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
532	! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
533	! , 2.tke_myj(i,1,j), akhs(i,j), akhs(i,j)ZSL_diag
534	! write(6,"(a)") &
535	! '{turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
536	! do k=kts,kte/2
537	! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
538	! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
539	! '{turb5 ', k, .01pmid(i,k,j),.01pint(i,k,j), T(i,k,j)-273.15 &
540	! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
541	! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
542	! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
543	! enddo
544	!
545	! ELSEIF(PRINT_DIAG==2)THEN
546	!
547	! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
548	! '}turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
549	! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
550	! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
551	! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
552	! '}turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
553	! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
554	! , 2.tke_myj(i,1,j), akhs(i,j), akhs(i,j)ZSL_diag
555	! write(6,"(a)") &
556	! '}turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
557	! do k=kts,kte/2
558	! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
559	! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
560	! '}turb5 ', k, .01pmid(i,k,j),.01pint(i,k,j), T(i,k,j)-273.15 &
561	! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
562	! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
563	! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
564	! enddo
565	! ENDIF
566	!*** End debugging
567	!
568	!----------------------------------------------------------------------
569	ENDDO
570	!----------------------------------------------------------------------
571	DO I=ITS,ITE
572	!
573	!*** FILL 1-D VERTICAL ARRAYS
574	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
575	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
576	!
577	DO K=KTS,KTE-1
578	AKMK(K)=AKM(I,K,J)
579	AKMK(K)=AKMK(K)(RHOK(K,I)+RHOK(K+1,I))0.5
580	ENDDO
581	!
582	LLOW=LOWLYR(I,J)
583	AKMS_DENS=AKMS(I,J)*RHOK(KTE+1-LLOW,I)
584	!
585	DO K=KTE,KTS,-1
586	KFLIP=KTE+1-K
587	UK(K)=U(I,KFLIP,J)
588	VK(K)=V(I,KFLIP,J)
589	ZHK(K)=ZINT(I,K,J)
590	ENDDO
591	ZHK(KTE+1)=ZINT(I,KTE+1,J)
592	!
593	!----------------------------------------------------------------------
594	!*** CARRY OUT THE VERTICAL DIFFUSION OF
595	!*** VELOCITY COMPONENTS
596	!----------------------------------------------------------------------
597	!
598	CALL VDIFV(LMH,DTDIF,UZ0(I,J),VZ0(I,J) &
599	& ,AKMS_DENS,UK,VK,AKMK,ZHK,RHOK(KTS,I) &
600	& ,IDS,IDE,JDS,JDE,KDS,KDE &
601	& ,IMS,IME,JMS,JME,KMS,KME &
602	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
603	!
604	!----------------------------------------------------------------------
605	!***
606	!*** COMPUTE PRIMARY VARIABLE TENDENCIES
607	!***
608	DO K=KTS,KTE
609	KFLIP=KTE+1-K
610	DUDT=(UK(KFLIP)-U(I,K,J))*RDTTURBL
611	DVDT=(VK(KFLIP)-V(I,K,J))*RDTTURBL
612	RUBLTEN(I,K,J)=DUDT
613	RVBLTEN(I,K,J)=DVDT
614	ENDDO
615	!
616	ENDDO
617	!----------------------------------------------------------------------
618	!
619	ENDDO main_integration
620	!
621	!----------------------------------------------------------------------
622	!
623	END SUBROUTINE MYJPBL
624	!
625	!----------------------------------------------------------------------
626	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
627	!----------------------------------------------------------------------
628	SUBROUTINE MIXLEN &
629	!----------------------------------------------------------------------
630	! ******************************************************************
631	! * *
632	! * LEVEL 2.5 MIXING LENGTH *
633	! * *
634	! ******************************************************************
635	!
636	&(LMH,U,V,T,THE,Q,CWM,Q2,Z,GM,GH,EL,PBLH,LPBL,LMXL,CT &
637	&,IDS,IDE,JDS,JDE,KDS,KDE &
638	&,IMS,IME,JMS,JME,KMS,KME &
639	&,ITS,ITE,JTS,JTE,KTS,KTE)
640	!----------------------------------------------------------------------
641	!
642	IMPLICIT NONE
643	!
644	!----------------------------------------------------------------------
645	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
646	& ,IMS,IME,JMS,JME,KMS,KME &
647	& ,ITS,ITE,JTS,JTE,KTS,KTE
648	!
649	INTEGER,INTENT(IN) :: LMH
650	!
651	INTEGER,INTENT(OUT) :: LMXL,LPBL
652	!
653	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: CWM,Q,Q2,T,THE,U,V
654	!
655	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
656	!
657	REAL,INTENT(OUT) :: PBLH
658	!
659	REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: EL,GH,GM
660	!
661	REAL,INTENT(INOUT) :: CT
662	!----------------------------------------------------------------------
663	!***
664	!*** LOCAL VARIABLES
665	!***
666	INTEGER :: K,LPBLM
667	!
668	REAL :: A,ADEN,B,BDEN,AUBR,BUBR,BLMX,EL0,ELOQ2X,GHL,GML &
669	& ,QOL2ST,QOL2UN,QDZL,RDZ,SQ,SREL,SZQ,TEM,THM,VKRMZ
670	!
671	REAL,DIMENSION(KTS:KTE) :: Q1
672	!
673	REAL,DIMENSION(KTS:KTE-1) :: DTH,ELM,REL
674	!
675	!----------------------------------------------------------------------
676	!**********************************************************************
677	!--------------FIND THE HEIGHT OF THE PBL-------------------------------
678	LPBL=LMH
679	!
680	DO K=LMH-1,1,-1
681	IF(Q2(K)<=EPSQ2*FH)THEN
682	LPBL=K
683	GO TO 110
684	ENDIF
685	ENDDO
686	!
687	LPBL=1
688	!
689	!--------------THE HEIGHT OF THE PBL------------------------------------
690	!
691	110 PBLH=Z(LPBL)-Z(LMH+1)
692	!
693	!-----------------------------------------------------------------------
694	DO K=KTS,LMH
695	Q1(K)=0.
696	ENDDO
697	!
698	DO K=1,LMH-1
699	DTH(K)=THE(K)-THE(K+1)
700	ENDDO
701	!
702	DO K=LMH-2,1,-1
703	IF(DTH(K)>0..AND.DTH(K+1)<=0.)THEN
704	DTH(K)=DTH(K)+CT
705	EXIT
706	ENDIF
707	ENDDO
708	!
709	CT=0.
710	!----------------------------------------------------------------------
711	DO K=KTS,LMH-1
712	RDZ=2./(Z(K)-Z(K+2))
713	GML=((U(K)-U(K+1))2+(V(K)-V(K+1))2)RDZRDZ
714	GM(K)=MAX(GML,EPSGM)
715	!
716	TEM=(T(K)+T(K+1))*0.5
717	THM=(THE(K)+THE(K+1))*0.5
718	!
719	A=THM*P608
720	B=(ELOCP/TEM-1.-P608)*THM
721	!
722	GHL=(DTH(K)((Q(K)+Q(K+1)+CWM(K)+CWM(K+1))(0.5*P608)+1.) &
723	& +(Q(K)-Q(K+1)+CWM(K)-CWM(K+1))*A &
724	& +(CWM(K)-CWM(K+1))B)RDZ
725	!
726	IF(ABS(GHL)<=EPSGH)GHL=EPSGH
727	GH(K)=GHL
728	ENDDO
729	!
730	!----------------------------------------------------------------------
731	!*** FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
732	!----------------------------------------------------------------------
733	!
734	LMXL=LMH
735	!
736	DO K=KTS,LMH-1
737	GML=GM(K)
738	GHL=GH(K)
739	!
740	IF(GHL>=EPSGH)THEN
741	IF(GML/GHL<=REQU)THEN
742	ELM(K)=EPSL
743	LMXL=K
744	ELSE
745	AUBR=(AUBMGML+AUBHGHL)*GHL
746	BUBR= BUBMGML+BUBHGHL
747	QOL2ST=(-0.5BUBR+SQRT(BUBRBUBR0.25-AUBRCUBR))*RCUBR
748	ELOQ2X=1./QOL2ST
749	ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
750	ENDIF
751	ELSE
752	ADEN=(ADNMGML+ADNHGHL)*GHL
753	BDEN= BDNMGML+BDNHGHL
754	QOL2UN=-0.5BDEN+SQRT(BDENBDEN*0.25-ADEN)
755	ELOQ2X=1./(QOL2UN+EPSRU) ! repsr1/qol2un
756	ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
757	ENDIF
758	ENDDO
759	!
760	IF(ELM(LMH-1)==EPSL)LMXL=LMH
761	!
762	!----------------------------------------------------------------------
763	!*** THE HEIGHT OF THE MIXED LAYER
764	!----------------------------------------------------------------------
765	!
766	BLMX=Z(LMXL)-Z(LMH+1)
767	!
768	!----------------------------------------------------------------------
769	DO K=LPBL,LMH
770	Q1(K)=SQRT(Q2(K))
771	ENDDO
772	!----------------------------------------------------------------------
773	SZQ=0.
774	SQ =0.
775	!
776	DO K=KTS,LMH-1
777	QDZL=(Q1(K)+Q1(K+1))*(Z(K+1)-Z(K+2))
778	SZQ=(Z(K+1)+Z(K+2)-Z(LMH+1)-Z(LMH+1))*QDZL+SZQ
779	SQ=QDZL+SQ
780	ENDDO
781	!
782	!----------------------------------------------------------------------
783	!*** COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
784	!----------------------------------------------------------------------
785	!
786	EL0=MIN(ALPHSZQ0.5/SQ,EL0MAX)
787	EL0=MAX(EL0 ,EL0MIN)
788	!
789	!----------------------------------------------------------------------
790	!*** ABOVE THE PBL TOP
791	!----------------------------------------------------------------------
792	!
793	LPBLM=MAX(LPBL-1,1)
794	!
795	DO K=KTS,LPBLM
796	EL(K)=MIN((Z(K)-Z(K+2))*ELFC,ELM(K))
797	REL(K)=EL(K)/ELM(K)
798	ENDDO
799	!
800	!----------------------------------------------------------------------
801	!*** INSIDE THE PBL
802	!----------------------------------------------------------------------
803	!
804	IF(LPBL<LMH)THEN
805	DO K=LPBL,LMH-1
806	VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
807	EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))
808	REL(K)=EL(K)/ELM(K)
809	ENDDO
810	ENDIF
811	!
812	DO K=LPBL+1,LMH-2
813	SREL=MIN(((REL(K-1)+REL(K+1))0.5+REL(K))0.5,REL(K))
814	EL(K)=MAX(SREL*ELM(K),EPSL)
815	ENDDO
816	!
817	!----------------------------------------------------------------------
818	END SUBROUTINE MIXLEN
819	!----------------------------------------------------------------------
820	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
821	!----------------------------------------------------------------------
822	SUBROUTINE PRODQ2 &
823	!----------------------------------------------------------------------
824	! ******************************************************************
825	! * *
826	! * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
827	! * *
828	! ******************************************************************
829	!
830	&(LMH,DTTURBL,USTAR,GM,GH,EL,Q2 &
831	&,IDS,IDE,JDS,JDE,KDS,KDE &
832	&,IMS,IME,JMS,JME,KMS,KME &
833	&,ITS,ITE,JTS,JTE,KTS,KTE)
834	!----------------------------------------------------------------------
835	!
836	IMPLICIT NONE
837	!
838	!----------------------------------------------------------------------
839	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
840	& ,IMS,IME,JMS,JME,KMS,KME &
841	& ,ITS,ITE,JTS,JTE,KTS,KTE
842	!
843	INTEGER,INTENT(IN) :: LMH
844	!
845	REAL,INTENT(IN) :: DTTURBL,USTAR
846	!
847	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: GH,GM
848	REAL,DIMENSION(KTS:KTE-1),INTENT(INOUT) :: EL
849	!
850	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
851	!----------------------------------------------------------------------
852	!***
853	!*** LOCAL VARIABLES
854	!***
855	INTEGER :: K
856	!
857	REAL :: ADEN,AEQU,ANUM,ARHS,BDEN,BEQU,BNUM,BRHS,CDEN,CRHS &
858	& ,DLOQ1,ELOQ11,ELOQ12,ELOQ13,ELOQ21,ELOQ22,ELOQ31,ELOQ32 &
859	& ,ELOQ41,ELOQ42,ELOQ51,ELOQ52,ELOQN,EQOL2,GHL,GML &
860	& ,RDEN1,RDEN2,RHS2,RHSP1,RHSP2,RHST2
861	!
862	!----------------------------------------------------------------------
863	!**********************************************************************
864	!----------------------------------------------------------------------
865	!
866	main_integration: DO K=1,LMH-1
867	GML=GM(K)
868	GHL=GH(K)
869	!
870	!----------------------------------------------------------------------
871	!*** COEFFICIENTS OF THE EQUILIBRIUM EQUATION
872	!----------------------------------------------------------------------
873	!
874	AEQU=(AEQMGML+AEQHGHL)*GHL
875	BEQU= BEQMGML+BEQHGHL
876	!
877	!----------------------------------------------------------------------
878	!*** EQUILIBRIUM SOLUTION FOR L/Q
879	!----------------------------------------------------------------------
880	!
881	EQOL2=-0.5BEQU+SQRT(BEQUBEQU*0.25-AEQU)
882	!
883	!----------------------------------------------------------------------
884	!*** IS THERE PRODUCTION/DISSIPATION ?
885	!----------------------------------------------------------------------
886	!
887	IF((GML+GHL*GHL<=EPSTRB) &
888	& .OR.(GHL>=EPSGH.AND.GML/GHL<=REQU) &
889	& .OR.(EQOL2<=EPS2))THEN
890	!
891	!----------------------------------------------------------------------
892	!*** NO TURBULENCE
893	!----------------------------------------------------------------------
894	!
895	Q2(K)=EPSQ2
896	EL(K)=EPSL
897	!----------------------------------------------------------------------
898	!
899	ELSE
900	!
901	!----------------------------------------------------------------------
902	!*** TURBULENCE
903	!----------------------------------------------------------------------
904	!----------------------------------------------------------------------
905	!*** COEFFICIENTS OF THE TERMS IN THE NUMERATOR
906	!----------------------------------------------------------------------
907	!
908	ANUM=(ANMMGML+ANMHGHL)*GHL
909	BNUM= BNMMGML+BNMHGHL
910	!
911	!----------------------------------------------------------------------
912	!*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
913	!----------------------------------------------------------------------
914	!
915	ADEN=(ADNMGML+ADNHGHL)*GHL
916	BDEN= BDNMGML+BDNHGHL
917	CDEN= 1.
918	!
919	!----------------------------------------------------------------------
920	!*** COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.
921	!----------------------------------------------------------------------
922	!
923	ARHS=-(ANUMBDEN-BNUMADEN)*2.
924	BRHS=- ANUM*4.
925	CRHS=- BNUM*2.
926	!
927	!----------------------------------------------------------------------
928	!*** INITIAL VALUE OF L/Q
929	!----------------------------------------------------------------------
930	!
931	DLOQ1=EL(K)/SQRT(Q2(K))
932	!
933	!----------------------------------------------------------------------
934	!*** FIRST ITERATION FOR L/Q, RHS=0
935	!----------------------------------------------------------------------
936	!
937	ELOQ21=1./EQOL2
938	ELOQ11=SQRT(ELOQ21)
939	ELOQ31=ELOQ21*ELOQ11
940	ELOQ41=ELOQ21*ELOQ21
941	ELOQ51=ELOQ21*ELOQ31
942	!
943	!----------------------------------------------------------------------
944	!*** 1./DENOMINATOR
945	!----------------------------------------------------------------------
946	!
947	RDEN1=1./(ADENELOQ41+BDENELOQ21+CDEN)
948	!
949	!----------------------------------------------------------------------
950	!*** D(RHS)/D(L/Q)
951	!----------------------------------------------------------------------
952	!
953	RHSP1=(ARHSELOQ51+BRHSELOQ31+CRHSELOQ11)RDEN1*RDEN1
954	!
955	!----------------------------------------------------------------------
956	!*** FIRST-GUESS SOLUTION
957	!----------------------------------------------------------------------
958	!
959	ELOQ12=ELOQ11+(DLOQ1-ELOQ11)EXP(RHSP1DTTURBL)
960	ELOQ12=MAX(ELOQ12,EPS1)
961	!
962	!----------------------------------------------------------------------
963	!*** SECOND ITERATION FOR L/Q
964	!----------------------------------------------------------------------
965	!
966	ELOQ22=ELOQ12*ELOQ12
967	ELOQ32=ELOQ22*ELOQ12
968	ELOQ42=ELOQ22*ELOQ22
969	ELOQ52=ELOQ22*ELOQ32
970	!
971	!----------------------------------------------------------------------
972	!*** 1./DENOMINATOR
973	!----------------------------------------------------------------------
974	!
975	RDEN2=1./(ADENELOQ42+BDENELOQ22+CDEN)
976	RHS2 =-(ANUMELOQ42+BNUMELOQ22)*RDEN2+RB1
977	RHSP2= (ARHSELOQ52+BRHSELOQ32+CRHSELOQ12)RDEN2*RDEN2
978	RHST2=RHS2/RHSP2
979	!
980	!----------------------------------------------------------------------
981	!*** CORRECTED SOLUTION
982	!----------------------------------------------------------------------
983	!
984	ELOQ13=ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)EXP(RHSP2DTTURBL)
985	ELOQ13=AMAX1(ELOQ13,EPS1)
986	!
987	!----------------------------------------------------------------------
988	!*** TWO ITERATIONS IS ENOUGH IN MOST CASES ...
989	!----------------------------------------------------------------------
990	!
991	ELOQN=ELOQ13
992	!
993	IF(ELOQN>EPS1)THEN
994	Q2(K)=EL(K)EL(K)/(ELOQNELOQN)
995	Q2(K)=AMAX1(Q2(K),EPSQ2)
996	IF(Q2(K)==EPSQ2)THEN
997	EL(K)=EPSL
998	ENDIF
999	ELSE
1000	Q2(K)=EPSQ2
1001	EL(K)=EPSL
1002	ENDIF
1003	!
1004	!----------------------------------------------------------------------
1005	!*** END OF TURBULENT BRANCH
1006	!----------------------------------------------------------------------
1007	!
1008	ENDIF
1009	!----------------------------------------------------------------------
1010	!*** END OF PRODUCTION/DISSIPATION LOOP
1011	!----------------------------------------------------------------------
1012	!
1013	ENDDO main_integration
1014	!
1015	!----------------------------------------------------------------------
1016	!*** LOWER BOUNDARY CONDITION FOR Q2
1017	!----------------------------------------------------------------------
1018	!
1019	Q2(LMH)=AMAX1(B1*(2./3.)USTAR*USTAR,EPSQ2)
1020	!----------------------------------------------------------------------
1021	!
1022	END SUBROUTINE PRODQ2
1023	!
1024	!----------------------------------------------------------------------
1025	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1026	!----------------------------------------------------------------------
1027	SUBROUTINE DIFCOF &
1028	! ******************************************************************
1029	! * *
1030	! * LEVEL 2.5 DIFFUSION COEFFICIENTS *
1031	! * *
1032	! ******************************************************************
1033	&(LMH,LMXL,GM,GH,EL,T,Q2,Z,AKM,AKH &
1034	&,IDS,IDE,JDS,JDE,KDS,KDE &
1035	&,IMS,IME,JMS,JME,KMS,KME &
1036	&,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
1037	!----------------------------------------------------------------------
1038	!
1039	IMPLICIT NONE
1040	!
1041	!----------------------------------------------------------------------
1042	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1043	& ,IMS,IME,JMS,JME,KMS,KME &
1044	& ,ITS,ITE,JTS,JTE,KTS,KTE
1045	!
1046	INTEGER,INTENT(IN) :: LMH,LMXL
1047	!
1048	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: Q2,T
1049	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL,GH,GM
1050	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1051	!
1052	REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: AKH,AKM
1053	!----------------------------------------------------------------------
1054	!***
1055	!*** LOCAL VARIABLES
1056	!***
1057	INTEGER :: K,KINV
1058	!
1059	REAL :: ADEN,AKMIN,BDEN,BESH,BESM,CDEN,D2T,ELL,ELOQ2,ELOQ4,ELQDZ &
1060	& ,ESH,ESM,GHL,GML,Q1L,RDEN,RDZ
1061	!
1062	!*** Begin debugging
1063	INTEGER,INTENT(IN) :: PRINT_DIAG
1064	! REAL :: D2Tmin
1065	!*** End debugging
1066	!
1067	!----------------------------------------------------------------------
1068	!**********************************************************************
1069	!----------------------------------------------------------------------
1070	!
1071	DO K=1,LMH-1
1072	ELL=EL(K)
1073	!
1074	ELOQ2=ELL*ELL/Q2(K)
1075	ELOQ4=ELOQ2*ELOQ2
1076	!
1077	GML=GM(K)
1078	GHL=GH(K)
1079	!
1080	!----------------------------------------------------------------------
1081	!*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
1082	!----------------------------------------------------------------------
1083	!
1084	ADEN=(ADNMGML+ADNHGHL)*GHL
1085	BDEN= BDNMGML+BDNHGHL
1086	CDEN= 1.
1087	!
1088	!----------------------------------------------------------------------
1089	!*** COEFFICIENTS FOR THE SM DETERMINANT
1090	!----------------------------------------------------------------------
1091	!
1092	BESM=BSMH*GHL
1093	!
1094	!----------------------------------------------------------------------
1095	!*** COEFFICIENTS FOR THE SH DETERMINANT
1096	!----------------------------------------------------------------------
1097	!
1098	BESH=BSHMGML+BSHHGHL
1099	!
1100	!----------------------------------------------------------------------
1101	!*** 1./DENOMINATOR
1102	!----------------------------------------------------------------------
1103	!
1104	RDEN=1./(ADENELOQ4+BDENELOQ2+CDEN)
1105	!
1106	!----------------------------------------------------------------------
1107	!*** SM AND SH
1108	!----------------------------------------------------------------------
1109	!
1110	ESM=(BESMELOQ2+CESM)RDEN
1111	ESH=(BESHELOQ2+CESH)RDEN
1112	!
1113	!----------------------------------------------------------------------
1114	!*** DIFFUSION COEFFICIENTS
1115	!----------------------------------------------------------------------
1116	!
1117	RDZ=2./(Z(K)-Z(K+2))
1118	Q1L=SQRT(Q2(K))
1119	ELQDZ=ELLQ1LRDZ
1120	AKM(K)=ELQDZ*ESM
1121	AKH(K)=ELQDZ*ESH
1122	!----------------------------------------------------------------------
1123	ENDDO
1124	!----------------------------------------------------------------------
1125	!
1126	!----------------------------------------------------------------------
1127	!*** INVERSIONS
1128	!----------------------------------------------------------------------
1129	!
1130	! IF(LMXL==LMH)THEN
1131	! KINV=LMH
1132	! D2Tmin=0.
1133	!
1134	! DO K=LMH/2,LMH-1
1135	! D2T=T(K-1)-2.*T(K)+T(K+1)
1136	! IF(D2T<D2Tmin)THEN
1137	! D2Tmin=D2T
1138	! IF(D2T<0)KINV=K
1139	! ENDIF
1140	! ENDDO
1141	!
1142	! IF(KINV<LMH)THEN
1143	! DO K=KINV-1,LMH-1
1144	! RDZ=2./(Z(K)-Z(K+2))
1145	! AKMIN=0.5*RDZ
1146	! AKM(K)=MAX(AKM(K),AKMIN)
1147	! AKH(K)=MAX(AKH(K),AKMIN)
1148	! ENDDO
1149	!
1150	!*** Begin debugging
1151	! IF(PRINT_DIAG>0)THEN
1152	! write(6,"(a,3i3)") '{turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
1153	! write(6,"(a,3i3)") '}turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
1154	! IF(PRINT_DIAG==1)THEN
1155	! write(6,"(a)") &
1156	! '{turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
1157	! ELSE
1158	! write(6,"(a)") &
1159	! '}turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
1160	! ENDIF
1161	! DO K=LMH-1,KINV-1,-1
1162	! D2T=T(K-1)-2.*T(K)+T(K+1)
1163	! RDZ=2./(Z(K)-Z(K+2))
1164	! AKMIN=0.5*RDZ
1165	! IF(PRINT_DIAG==1)THEN
1166	! write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '{turb3 ' &
1167	! ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
1168	! ELSE
1169	! write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '}turb3 ' &
1170	! ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
1171	! ENDIF
1172	! ENDDO
1173	! ENDIF !- IF (print_diag > 0) THEN
1174	! ENDIF !- IF(KINV<LMH)THEN
1175	!*** End debugging
1176	!
1177	! ENDIF
1178	!----------------------------------------------------------------------
1179	!
1180	END SUBROUTINE DIFCOF
1181	!
1182	!----------------------------------------------------------------------
1183	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1184	!----------------------------------------------------------------------
1185	SUBROUTINE VDIFQ &
1186	! ******************************************************************
1187	! * *
1188	! * VERTICAL DIFFUSION OF Q2 (TKE) *
1189	! * *
1190	! ******************************************************************
1191	&(LMH,DTDIF,Q2,EL,Z &
1192	&,IDS,IDE,JDS,JDE,KDS,KDE &
1193	&,IMS,IME,JMS,JME,KMS,KME &
1194	&,ITS,ITE,JTS,JTE,KTS,KTE)
1195	!----------------------------------------------------------------------
1196	!
1197	IMPLICIT NONE
1198	!
1199	!----------------------------------------------------------------------
1200	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1201	& ,IMS,IME,JMS,JME,KMS,KME &
1202	& ,ITS,ITE,JTS,JTE,KTS,KTE
1203	!
1204	INTEGER,INTENT(IN) :: LMH
1205	!
1206	REAL,INTENT(IN) :: DTDIF
1207	!
1208	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL
1209	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1210	!
1211	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
1212	!----------------------------------------------------------------------
1213	!***
1214	!*** LOCAL VARIABLES
1215	!***
1216	INTEGER :: K
1217	!
1218	REAL :: ADEN,AKQS,BDEN,BESH,BESM,CDEN,CF,DTOZS,ELL,ELOQ2,ELOQ4 &
1219	& ,ELQDZ,ESH,ESM,ESQHF,GHL,GML,Q1L,RDEN,RDZ
1220	!
1221	REAL,DIMENSION(KTS:KTE-2) :: AKQ,CM,CR,DTOZ,RSQ2
1222	!----------------------------------------------------------------------
1223	!**********************************************************************
1224	!----------------------------------------------------------------------
1225	!***
1226	!*** VERTICAL TURBULENT DIFFUSION
1227	!***
1228	!----------------------------------------------------------------------
1229	ESQHF=0.5*ESQ
1230	!
1231	DO K=KTS,LMH-2
1232	DTOZ(K)=(DTDIF+DTDIF)/(Z(K)-Z(K+2))
1233	AKQ(K)=SQRT((Q2(K)+Q2(K+1))0.5)(EL(K)+EL(K+1))*ESQHF &
1234	& /(Z(K+1)-Z(K+2))
1235	CR(K)=-DTOZ(K)*AKQ(K)
1236	ENDDO
1237	!
1238	CM(1)=DTOZ(1)*AKQ(1)+1.
1239	RSQ2(1)=Q2(1)
1240	!
1241	DO K=KTS+1,LMH-2
1242	CF=-DTOZ(K)*AKQ(K-1)/CM(K-1)
1243	CM(K)=-CR(K-1)CF+(AKQ(K-1)+AKQ(K))DTOZ(K)+1.
1244	RSQ2(K)=-RSQ2(K-1)*CF+Q2(K)
1245	ENDDO
1246	!
1247	DTOZS=(DTDIF+DTDIF)/(Z(LMH-1)-Z(LMH+1))
1248	AKQS=SQRT((Q2(LMH-1)+Q2(LMH))0.5)(EL(LMH-1)+ELZ0)*ESQHF &
1249	& /(Z(LMH)-Z(LMH+1))
1250	!
1251	CF=-DTOZS*AKQ(LMH-2)/CM(LMH-2)
1252	!
1253	Q2(LMH-1)=(DTOZSAKQSQ2(LMH)-RSQ2(LMH-2)*CF+Q2(LMH-1)) &
1254	& /((AKQ(LMH-2)+AKQS)DTOZS-CR(LMH-2)CF+1.)
1255	!
1256	DO K=LMH-2,KTS,-1
1257	Q2(K)=(-CR(K)*Q2(K+1)+RSQ2(K))/CM(K)
1258	ENDDO
1259	!----------------------------------------------------------------------
1260	!
1261	END SUBROUTINE VDIFQ
1262	!
1263	!----------------------------------------------------------------------
1264	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1265	!---------------------------------------------------------------------
1266	SUBROUTINE VDIFH(DTDIF,LMH,THZ0,QZ0,RKHS,CHKLOWQ,CT &
1267	& ,THE,Q,CWM,RKH,Z,RHO &
1268	& ,IDS,IDE,JDS,JDE,KDS,KDE &
1269	& ,IMS,IME,JMS,JME,KMS,KME &
1270	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1271	! ***************************************************************
1272	! * *
1273	! * VERTICAL DIFFUSION OF MASS VARIABLES *
1274	! * *
1275	! ***************************************************************
1276	!---------------------------------------------------------------------
1277	!
1278	IMPLICIT NONE
1279	!
1280	!---------------------------------------------------------------------
1281	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1282	& ,IMS,IME,JMS,JME,KMS,KME &
1283	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1284	!
1285	INTEGER,INTENT(IN) :: LMH
1286	!
1287	REAL,INTENT(IN) :: CHKLOWQ,CT,DTDIF,QZ0,RKHS,THZ0
1288	!
1289	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKH
1290	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1291	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1292	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: CWM,Q,THE
1293	!
1294	!----------------------------------------------------------------------
1295	!***
1296	!*** LOCAL VARIABLES
1297	!***
1298	INTEGER :: K
1299	!
1300	REAL :: CF,CMB,CMCB,CMQB,CMTB,CTHF,DTOZL,DTOZS &
1301	& ,RCML,RKHH,RKQS,RSCB,RSQB,RSTB
1302	!
1303	REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RKCT,RSC,RSQ,RST
1304	!
1305	!----------------------------------------------------------------------
1306	!**********************************************************************
1307	!----------------------------------------------------------------------
1308	CTHF=0.5*CT
1309	!
1310	DO K=KTS,LMH-1
1311	DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1312	CR(K)=-DTOZ(K)*RKH(K)
1313	RKCT(K)=RKH(K)(Z(K)-Z(K+2))CTHF
1314	ENDDO
1315	!
1316	CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHO(KTS)
1317	!----------------------------------------------------------------------
1318	RST(KTS)=-RKCT(KTS)*DTOZ(KTS) &
1319	& +THE(KTS)*RHO(KTS)
1320	RSQ(KTS)=Q(KTS) *RHO(KTS)
1321	RSC(KTS)=CWM(KTS)*RHO(KTS)
1322	!----------------------------------------------------------------------
1323	DO K=KTS+1,LMH-1
1324	DTOZL=DTOZ(K)
1325	CF=-DTOZL*RKH(K-1)/CM(K-1)
1326	CM(K)=-CR(K-1)CF+(RKH(K-1)+RKH(K))DTOZL+RHO(K)
1327	RST(K)=-RST(K-1)CF+(RKCT(K-1)-RKCT(K))DTOZL+THE(K)*RHO(K)
1328	RSQ(K)=-RSQ(K-1)CF+Q(K) RHO(K)
1329	RSC(K)=-RSC(K-1)CF+CWM(K)RHO(K)
1330	ENDDO
1331	!
1332	DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1333	RKHH=RKH(LMH-1)
1334	!
1335	CF=-DTOZS*RKHH/CM(LMH-1)
1336	RKQS=RKHS*CHKLOWQ
1337	!
1338	CMB=CR(LMH-1)*CF
1339	CMTB=-CMB+(RKHH+RKHS)*DTOZS+RHO(LMH)
1340	CMQB=-CMB+(RKHH+RKQS)*DTOZS+RHO(LMH)
1341	CMCB=-CMB+(RKHH )*DTOZS+RHO(LMH)
1342	!
1343	RSTB=-RST(LMH-1)CF+RKCT(LMH-1)DTOZS+THE(LMH)*RHO(LMH)
1344	RSQB=-RSQ(LMH-1)CF+Q(LMH) RHO(LMH)
1345	RSCB=-RSC(LMH-1)CF+CWM(LMH)RHO(LMH)
1346	!----------------------------------------------------------------------
1347	THE(LMH)=(DTOZSRKHSTHZ0+RSTB)/CMTB
1348	Q(LMH) =(DTOZSRKQSQZ0 +RSQB)/CMQB
1349	CWM(LMH)=( RSCB)/CMCB
1350	!----------------------------------------------------------------------
1351	DO K=LMH-1,KTS,-1
1352	RCML=1./CM(K)
1353	THE(K)=(-CR(K)THE(K+1)+RST(K))RCML
1354	Q(K) =(-CR(K)* Q(K+1)+RSQ(K))*RCML
1355	CWM(K)=(-CR(K)CWM(K+1)+RSC(K))RCML
1356	ENDDO
1357	!----------------------------------------------------------------------
1358	!
1359	END SUBROUTINE VDIFH
1360	!
1361	!----------------------------------------------------------------------
1362	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1363	!---------------------------------------------------------------------
1364	SUBROUTINE VDIFX(DTDIF,LMH,RKHS,CT &
1365	,DUST1,DUST2,RKH,Z,RHO &
1366	,IDS,IDE,JDS,JDE,KDS,KDE &
1367	,IMS,IME,JMS,JME,KMS,KME &
1368	,ITS,ITE,JTS,JTE,KTS,KTE)
1369	! ***************************************************************
1370	! * *
1371	! * VERTICAL DIFFUSION OF MASS VARIABLES *
1372	! * *
1373	! ***************************************************************
1374	!---------------------------------------------------------------------
1375	!
1376	IMPLICIT NONE
1377	!
1378	!---------------------------------------------------------------------
1379	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1380	,IMS,IME,JMS,JME,KMS,KME &
1381	,ITS,ITE,JTS,JTE,KTS,KTE
1382	!
1383	INTEGER,INTENT(IN) :: LMH
1384	!
1385	REAL,INTENT(IN) :: CT,DTDIF,RKHS
1386	!
1387	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKH
1388	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1389	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1390	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: DUST1,DUST2
1391	!
1392	!----------------------------------------------------------------------
1393	!***
1394	!*** LOCAL VARIABLES
1395	!***
1396	INTEGER :: K
1397	!
1398	REAL :: CF,CMB,CMDB,CTHF,DTOZL,DTOZS &
1399	,RCML,RKHH,RSD1B,RSD2B
1400	!
1401	REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RKCT,RSD1,RSD2
1402	!
1403	!----------------------------------------------------------------------
1404	!**********************************************************************
1405	!----------------------------------------------------------------------
1406	CTHF=0.5*CT
1407	!
1408	DO K=KTS,LMH-1
1409	DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1410	CR(K)=-DTOZ(K)*RKH(K)
1411	RKCT(K)=RKH(K)(Z(K)-Z(K+2))CTHF
1412	ENDDO
1413	!
1414	CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHO(KTS)
1415	!----------------------------------------------------------------------
1416	RSD1(KTS)=DUST1(KTS)*RHO(KTS)
1417	RSD2(KTS)=DUST2(KTS)*RHO(KTS)
1418	!----------------------------------------------------------------------
1419	DO K=KTS+1,LMH-1
1420	DTOZL=DTOZ(K)
1421	CF=-DTOZL*RKH(K-1)/CM(K-1)
1422	CM(K)=-CR(K-1)CF+(RKH(K-1)+RKH(K))DTOZL+RHO(K)
1423	RSD1(K)=-RSD1(K-1)CF+DUST1(K)RHO(K)
1424	RSD2(K)=-RSD2(K-1)CF+DUST2(K)RHO(K)
1425	ENDDO
1426	!
1427	DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1428	RKHH=RKH(LMH-1)
1429	!
1430	CF=-DTOZS*RKHH/CM(LMH-1)
1431	!
1432	CMB=CR(LMH-1)*CF
1433	CMDB=-CMB+RKHH*DTOZS+RHO(LMH)
1434	!
1435	RSD1B=-RSD1(LMH-1)CF+DUST1(LMH)RHO(LMH)
1436	RSD2B=-RSD2(LMH-1)CF+DUST2(LMH)RHO(LMH)
1437	!----------------------------------------------------------------------
1438	DUST1(LMH)=RSD1B/CMDB
1439	DUST2(LMH)=RSD2B/CMDB
1440	!----------------------------------------------------------------------
1441	DO K=LMH-1,KTS,-1
1442	RCML=1./CM(K)
1443	DUST1(K)=(-CR(K)DUST1(K+1)+RSD1(K))RCML
1444	DUST2(K)=(-CR(K)DUST2(K+1)+RSD2(K))RCML
1445	ENDDO
1446	!----------------------------------------------------------------------
1447	!
1448	END SUBROUTINE VDIFX
1449
1450	!---------------------------------------------------------------------
1451	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1452	!---------------------------------------------------------------------
1453	SUBROUTINE VDIFV(LMH,DTDIF,UZ0,VZ0,RKMS,U,V,RKM,Z,RHO &
1454	& ,IDS,IDE,JDS,JDE,KDS,KDE &
1455	& ,IMS,IME,JMS,JME,KMS,KME &
1456	,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1457	! ***************************************************************
1458	! * *
1459	! * VERTICAL DIFFUSION OF VELOCITY COMPONENTS *
1460	! * *
1461	! ***************************************************************
1462	!---------------------------------------------------------------------
1463	!
1464	IMPLICIT NONE
1465	!
1466	!---------------------------------------------------------------------
1467	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1468	& ,IMS,IME,JMS,JME,KMS,KME &
1469	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1470	!
1471	INTEGER,INTENT(IN) :: LMH
1472	!
1473	REAL,INTENT(IN) :: RKMS,DTDIF,UZ0,VZ0
1474	!
1475	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKM
1476	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1477	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1478	!
1479	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: U,V
1480	!----------------------------------------------------------------------
1481	!***
1482	!*** LOCAL VARIABLES
1483	!***
1484	INTEGER :: K
1485	!
1486	REAL :: CF,DTOZAK,DTOZL,DTOZS,RCML,RCMVB,RHOK,RKMH
1487	!
1488	REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RSU,RSV
1489	!----------------------------------------------------------------------
1490	!**********************************************************************
1491	!----------------------------------------------------------------------
1492	DO K=1,LMH-1
1493	DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1494	CR(K)=-DTOZ(K)*RKM(K)
1495	ENDDO
1496	!
1497	RHOK=RHO(1)
1498	CM(1)=DTOZ(1)*RKM(1)+RHOK
1499	RSU(1)=U(1)*RHOK
1500	RSV(1)=V(1)*RHOK
1501	!----------------------------------------------------------------------
1502	DO K=2,LMH-1
1503	DTOZL=DTOZ(K)
1504	CF=-DTOZL*RKM(K-1)/CM(K-1)
1505	RHOK=RHO(K)
1506	CM(K)=-CR(K-1)CF+(RKM(K-1)+RKM(K))DTOZL+RHOK
1507	RSU(K)=-RSU(K-1)CF+U(K)RHOK
1508	RSV(K)=-RSV(K-1)CF+V(K)RHOK
1509	ENDDO
1510	!----------------------------------------------------------------------
1511	DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1512	RKMH=RKM(LMH-1)
1513	!
1514	CF=-DTOZS*RKMH/CM(LMH-1)
1515	RHOK=RHO(LMH)
1516	RCMVB=1./((RKMH+RKMS)DTOZS-CR(LMH-1)CF+RHOK)
1517	DTOZAK=DTOZS*RKMS
1518	!----------------------------------------------------------------------
1519	U(LMH)=(DTOZAKUZ0-RSU(LMH-1)CF+U(LMH)RHOK)RCMVB
1520	V(LMH)=(DTOZAKVZ0-RSV(LMH-1)CF+V(LMH)RHOK)RCMVB
1521	!----------------------------------------------------------------------
1522	DO K=LMH-1,1,-1
1523	RCML=1./CM(K)
1524	U(K)=(-CR(K)U(K+1)+RSU(K))RCML
1525	V(K)=(-CR(K)V(K+1)+RSV(K))RCML
1526	ENDDO
1527	!----------------------------------------------------------------------
1528	!
1529	END SUBROUTINE VDIFV
1530	!
1531	!-----------------------------------------------------------------------
1532	!
1533	!=======================================================================
1534	SUBROUTINE MYJPBLINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1535	& TKE_MYJ,EXCH_H,RESTART,ALLOWED_TO_READ, &
1536	& IDS,IDE,JDS,JDE,KDS,KDE, &
1537	& IMS,IME,JMS,JME,KMS,KME, &
1538	& ITS,ITE,JTS,JTE,KTS,KTE )
1539	!-----------------------------------------------------------------------
1540	IMPLICIT NONE
1541	!-----------------------------------------------------------------------
1542	LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1543	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1544	& IMS,IME,JMS,JME,KMS,KME, &
1545	& ITS,ITE,JTS,JTE,KTS,KTE
1546
1547	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
1548	& RUBLTEN, &
1549	& RVBLTEN, &
1550	& RTHBLTEN, &
1551	& RQVBLTEN, &
1552	& TKE_MYJ
1553	INTEGER :: I,J,K,ITF,JTF,KTF
1554	!-----------------------------------------------------------------------
1555	!-----------------------------------------------------------------------
1556
1557	JTF=MIN0(JTE,JDE-1)
1558	KTF=MIN0(KTE,KDE-1)
1559	ITF=MIN0(ITE,IDE-1)
1560
1561	IF(.NOT.RESTART)THEN
1562	DO J=JTS,JTF
1563	DO K=KTS,KTF
1564	DO I=ITS,ITF
1565	TKE_MYJ(I,K,J)=EPSQ2
1566	RUBLTEN(I,K,J)=0.
1567	RVBLTEN(I,K,J)=0.
1568	RTHBLTEN(I,K,J)=0.
1569	RQVBLTEN(I,K,J)=0.
1570	EXCH_H(I,K,J)=0.
1571	ENDDO
1572	ENDDO
1573	ENDDO
1574	ENDIF
1575
1576	END SUBROUTINE MYJPBLINIT
1577	!-----------------------------------------------------------------------
1578	!
1579	END MODULE MODULE_BL_MYJPBL
1580	!
1581	!-----------------------------------------------------------------------

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

Download in other formats:

Original Format