Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_myjpbl.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 63.1 KB

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

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

Download in other formats:

Original Format