Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_qnsepbl.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: 50.0 KB

Line
1	!-----------------------------------------------------------------------
2	!
3	MODULE MODULE_BL_QNSEPBL
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	!
31	!-----------------------------------------------------------------------
32	!*** QNSE MODEL CONSTANTS
33	!-----------------------------------------------------------------------
34	!
35	REAL,PARAMETER :: EPSQ2L=0.01
36	REAL,PARAMETER :: C0=0.55,CEPS=C0**3,BLCKDR=0.0063,CN=0.75 &
37	& ,AM1=8.0,AM2=2.3,AM3=35.0,AH1=1.4,AH2=-0.01 &
38	& ,AH3=1.29,AH4=2.44,AH5=19.8 &
39	& ,ARIMIN=0.127,BM1=2.88,BM2=16.0,BH1=3.6,BH2=16.0 &
40	& ,BH3=720.0,EPSKM=1.E-3
41	!
42	REAL,PARAMETER :: CAPA=R_D/CP
43	REAL,PARAMETER :: RLIVWV=XLS/XLV,ELOCP=2.72E6/CP
44	REAL,PARAMETER :: EPS1=1.E-12,EPS2=0.
45	REAL,PARAMETER :: EPSL=0.32,EPSRU=1.E-7,EPSRS=1.E-7 &
46	& ,EPSTRB=1.E-24
47	REAL,PARAMETER :: EPSA=1.E-8,EPSIT=1.E-4,EPSU2=1.E-4,EPSUST=0.07 &
48	& ,FH=1.01
49	REAL,PARAMETER :: ALPH=0.30,BETA=1./273.,EL0MAX=1000.,EL0MIN=1. &
50	& ,ELFC=0.23*0.5,GAM1=0.2222222222222222222 &
51	& ,PRT=1.
52	REAL,PARAMETER :: A1=0.659888514560862645 &
53	& ,A2x=0.6574209922667784586 &
54	& ,B1=11.87799326209552761 &
55	& ,B2=7.226971804046074028 &
56	& ,C1=0.000830955950095854396
57	REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
58	REAL,PARAMETER :: ELZ0=0.,ESQ=5.0,EXCM=0.001 &
59	& ,FHNEU=0.8,GLKBR=10.,GLKBS=30. &
60	& ,QVISC=2.1E-5,RFC=0.191,RIC=0.505,SMALL=0.35 &
61	& ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5 &
62	& ,USTC=0.7,USTR=0.225,VISC=1.5E-5 &
63	& ,WOLD=0.15,WWST=1.2,ZTMAX=1.,ZTFC=1.,ZTMIN=-5.
64	!
65	REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
66	!
67	REAL,PARAMETER :: BTG=BETAG,CZIV=SMALLGLKBS &
68	! & ,EP_1=R_V/R_D-1.,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
69	& ,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
70	& ,RB1=1./B1,RTVISC=1./TVISC,RVISC=1./VISC &
71	& ,ZQRZT=SQSC/SQPR
72	!
73	REAL,PARAMETER :: ADNH= 9.A1A2xA2x(12.A1+3.B2)BTGBTG &
74	& ,ADNM=18.A1A1A2x(B2-3.A2x)BTG &
75	& ,ANMH=-9.A1A2xA2xBTG*BTG &
76	& ,ANMM=-3.A1A2x(3.A2x+3.B2C1+18.A1C1-B2) &
77	& *BTG &
78	& ,BDNH= 3.A2x(7.A1+B2)BTG &
79	& ,BDNM= 6.A1A1 &
80	& ,BEQH= A2xB1BTG+3.A2x(7.A1+B2)BTG &
81	& ,BEQM=-A1B1(1.-3.C1)+6.A1*A1 &
82	& ,BNMH=-A2x*BTG &
83	& ,BNMM=A1(1.-3.C1) &
84	& ,BSHH=9.A1A2xA2xBTG &
85	& ,BSHM=18.A1A1A2xC1 &
86	& ,BSMH=-3.A1A2x(3.A2x+3.B2C1+12.A1C1-B2) &
87	& *BTG &
88	& ,CESH=A2x &
89	& ,CESM=A1(1.-3.C1) &
90	& ,CNV=EP_1*G/BTG &
91	& ,ELFCS=VKARMAN*BTG &
92	& ,FZQ1=RTVISCQVISCZQRZT &
93	& ,FZQ2=RTVISCQVISCZQRZT &
94	& ,FZT1=RVISC TVISCSQPR &
95	& ,FZT2=CZIVGRRSTVISC*SQPR &
96	& ,FZU1=CZIV*VISC &
97	& ,PIHF=0.5*PI &
98	& ,RFAC=RIC/(FHNEURFCRFC) &
99	& ,RQVISC=1./QVISC &
100	& ,RRIC=1./RIC &
101	& ,USTFC=0.018/G &
102	& ,WNEW=1.-WOLD &
103	& ,WWST2=WWST*WWST
104	!
105	!-----------------------------------------------------------------------
106	!* FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)2
107	!-----------------------------------------------------------------------
108	!
109	REAL,PARAMETER :: AEQH=9.A1A2xA2xB1BTGBTG &
110	& +9.A1A2xA2x(12.A1+3.B2)BTGBTG &
111	& ,AEQM=3.A1A2xB1(3.A2x+3.B2C1+18.A1*C1-B2)&
112	& BTG+18.A1A1A2x(B2-3.A2x)*BTG
113	!
114	!-----------------------------------------------------------------------
115	!*** FORBIDDEN TURBULENCE AREA
116	!-----------------------------------------------------------------------
117	!
118	REAL,PARAMETER :: REQU=-AEQH/AEQM &
119	& ,EPSGH=1.E-9,EPSGM=REQU*EPSGH
120	!
121	!-----------------------------------------------------------------------
122	!*** NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
123	!-----------------------------------------------------------------------
124	!
125	REAL,PARAMETER :: UBRYL=(18.REQUA1A1A2xB2C1*BTG &
126	& +9.A1A2xA2xB2BTGBTG) &
127	& /(REQU*ADNM+ADNH) &
128	& ,UBRY=(1.+EPSRS)UBRYL,UBRY3=3.UBRY
129	!
130	REAL,PARAMETER :: AUBH=27.A1A2xA2xB2BTGBTG-ADNH*UBRY3 &
131	& ,AUBM=54.A1A1A2xB2C1BTG -ADNM*UBRY3 &
132	& ,BUBH=(9.A1A2x+3.A2xB2)BTG-BDNHUBRY3 &
133	& ,BUBM=18.A1A1C1 -BDNMUBRY3 &
134	& ,CUBR=1. - UBRY3 &
135	& ,RCUBR=1./CUBR
136	!
137	!-----------------------------------------------------------------------
138	!
139	CONTAINS
140	!
141	!----------------------------------------------------------------------
142	SUBROUTINE QNSEPBL(DT,STEPBL,HT,DZ &
143	& ,PMID,PINT,TH,T,EXNER,QV,CWM,U,V,RHO &
144	& ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0,CORF &
145	& ,LOWLYR,XLAND,SICE,SNOW &
146	& ,TKE,EXCH_H,EXCH_M,USTAR,ZNT,EL_MYJ,PBLH,KPBL,CT &
147	& ,AKHS,AKMS,ELFLX &
148	& ,RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN &
149	& ,IDS,IDE,JDS,JDE,KDS,KDE &
150	& ,IMS,IME,JMS,JME,KMS,KME &
151	& ,ITS,ITE,JTS,JTE,KTS,KTE)
152	!----------------------------------------------------------------------
153	!
154	IMPLICIT NONE
155	!
156	!----------------------------------------------------------------------
157	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
158	& ,IMS,IME,JMS,JME,KMS,KME &
159	& ,ITS,ITE,JTS,JTE,KTS,KTE
160	!
161	INTEGER,INTENT(IN) :: STEPBL
162
163	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
164	!
165	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: KPBL
166	!
167	REAL,INTENT(IN) :: DT
168	!
169	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,SICE,SNOW &
170	& ,TSK,XLAND
171	!
172	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DZ &
173	& ,EXNER &
174	& ,PMID,PINT &
175	& ,QV,RHO &
176	& ,T,TH,U,V
177	!
178	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PBLH
179	!
180	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS
181	!
182	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
183	& ,INTENT(OUT) :: EL_MYJ &
184	& ,RQCBLTEN,RQVBLTEN &
185	& ,RTHBLTEN &
186	& ,RUBLTEN,RVBLTEN
187	!
188	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CT,QSFC,QZ0 &
189	& ,THZ0,USTAR &
190	& ,UZ0,VZ0,ZNT
191	!
192	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
193	& ,INTENT(INOUT) :: EXCH_H,EXCH_M,TKE
194	!
195	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHKLOWQ,ELFLX,CORF
196	!
197	!----------------------------------------------------------------------
198	!***
199	!*** LOCAL VARIABLES
200	!***
201	INTEGER :: I,J,K,KFLIP,LLOW,LMH,LMXL
202	!
203	INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: LPBL
204	!
205	REAL :: AKHS_DENS,AKMS_DENS,APEX,DCDT,DELTAZ,DQDT,DTDIF,DTDT &
206	& ,DTTURBL,DUDT,DVDT,EXNSFC,PSFC,PTOP,QFC1,QLOW,QOLD &
207	& ,RATIOMX,RDTTURBL,RG,RWMSK,SEAMASK,THNEW,THOLD,TX &
208	& ,ULOW,VLOW,WMSK
209	!
210	REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,TK,UK,VK
211	!
212	REAL,DIMENSION(KTS:KTE-1) :: AKHK,AKMK,EL,RI,GH,S2
213	!
214	REAL,DIMENSION(KTS:KTE+1) :: ZHK
215	!
216	REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
217	!
218	REAL,DIMENSION(KTS:KTE,ITS:ITE) :: RHOK
219	!
220	REAL,DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: APE,THE
221	!
222	REAL,DIMENSION(ITS:ITE,KTS:KTE-1,JTS:JTE) :: AKH,AKM
223	!
224	REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
225	!
226	!*** Begin debugging
227	REAL :: ZSL_DIAG
228	INTEGER :: IMD,JMD,PRINT_DIAG
229	!*** End debugging
230	!
231	!----------------------------------------------------------------------
232	!**********************************************************************
233	!----------------------------------------------------------------------
234	!
235	!*** Begin debugging
236	IMD=(IMS+IME)/2
237	JMD=(JMS+JME)/2
238	!*** End debugging
239	!
240	!*** MAKE PREPARATIONS
241	!
242	!----------------------------------------------------------------------
243	DTTURBL=DT*STEPBL
244	RDTTURBL=1./DTTURBL
245	DTDIF=DTTURBL
246	RG=1./G
247	!
248	DO J=JTS,JTE
249	DO K=KTS,KTE-1
250	DO I=ITS,ITE
251	AKM(I,K,J)=0.
252	ENDDO
253	ENDDO
254	ENDDO
255	!
256	DO J=JTS,JTE
257	DO K=KTS,KTE+1
258	DO I=ITS,ITE
259	ZINT(I,K,J)=0.
260	ENDDO
261	ENDDO
262	ENDDO
263	!
264	DO J=JTS,JTE
265	DO I=ITS,ITE
266	ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer
267	!
268	!!!!!!!!!
269	!!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
270	!!!!!!!!!
271	!!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer
272	!!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer
273	!
274	ENDDO
275	ENDDO
276	!
277	DO J=JTS,JTE
278	DO K=KTE,KTS,-1
279	KFLIP=KTE+1-K
280	DO I=ITS,ITE
281	ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
282	APEX=1./EXNER(I,K,J)
283	APE(I,K,J)=APEX
284	TX=T(I,K,J)
285	THE(I,K,J)=(CWM(I,K,J)(-ELOCP/TX)+1.)TH(I,K,J)
286	ENDDO
287	ENDDO
288	ENDDO
289	!
290	EL_MYJ(its:ite,:,jts:jte) = 0.
291	!
292	!----------------------------------------------------------------------
293	setup_integration: DO J=JTS,JTE
294	!----------------------------------------------------------------------
295	!
296	DO I=ITS,ITE
297	!
298	!*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
299	!
300	LMH=KTE-LOWLYR(I,J)+1
301	!
302	PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
303	PSFC=PINT(I,LOWLYR(I,J),J)
304	!
305	!*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
306	!
307	SEAMASK=XLAND(I,J)-1.
308	!
309	!*** FILL 1-D VERTICAL ARRAYS
310	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
311	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
312	!
313	DO K=KTE,KTS,-1
314	KFLIP=KTE+1-K
315	TK(K)=T(I,KFLIP,J)
316	THEK(K)=THE(I,KFLIP,J)
317	RATIOMX=QV(I,KFLIP,J)
318	QK(K)=RATIOMX/(1.+RATIOMX)
319	CWMK(K)=CWM(I,KFLIP,J)
320	PK(K)=PMID(I,KFLIP,J)
321	UK(K)=U(I,KFLIP,J)
322	VK(K)=V(I,KFLIP,J)
323	!
324	!*** TKE=0.5(q2) ==> q2=2.TKE
325	!
326	Q2K(K)=2.*TKE(I,KFLIP,J)
327	!
328	!*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
329	!
330	ZHK(K)=ZINT(I,K,J)
331	!
332	ENDDO
333	ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer
334	!
335	!*** Begin debugging
336	! IF(I==IMD.AND.J==JMD)THEN
337	! PRINT_DIAG=1
338	! ELSE
339	! PRINT_DIAG=0
340	! ENDIF
341	! IF(I==227.AND.J==363)PRINT_DIAG=2
342	!*** End debugging
343	!
344	!----------------------------------------------------------------------
345	!***
346	!*** FIND THE MIXING LENGTH
347	!***
348	CALL MIXLEN(LMH,UK,VK,TK,THEK,QK,CWMK &
349	& ,Q2K,ZHK,USTAR(I,J),CORF(I,J),S2,GH,RI,EL &
350	& ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J) &
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 QNSE EXCHANGE COEFFICIENTS
364	!***
365	CALL DIFCOF(LMH,EL,RI,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	!----------------------------------------------------------------------
371	!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
372	!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. COUNTING
373	!*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
374	!*** ARE DEFINED ON THE TOPS OF THE LAYERS KTS TO KTE-1.
375	!
376	DO K=KTS,KTE-1
377	KFLIP=KTE-K
378	AKH(I,K,J)=AKHK(K)
379	AKM(I,K,J)=AKMK(K)
380	DELTAZ=0.5*(ZHK(KFLIP)-ZHK(KFLIP+2))
381	EXCH_H(I,K+1,J)=AKHK(KFLIP)*DELTAZ
382	EXCH_M(I,K+1,J)=AKMK(KFLIP)*DELTAZ
383	ENDDO
384	!
385	!----------------------------------------------------------------------
386	!***
387	!*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
388	!*** THE TURBULENT KINETIC ENERGY
389	!***
390	!
391	CALL PRODQ2(LMH,DTTURBL,USTAR(I,J),S2,RI,Q2K,EL,ZHK,AKMK,AKHK &
392	& ,IDS,IDE,JDS,JDE,KDS,KDE &
393	& ,IMS,IME,JMS,JME,KMS,KME &
394	& ,ITS,ITE,JTS,JTE,KTS,KTE)
395	!
396	!----------------------------------------------------------------------
397	!***
398	!*** CARRY OUT THE VERTICAL DIFFUSION OF
399	!*** TURBULENT KINETIC ENERGY
400	!***
401	!
402	CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK &
403	& ,IDS,IDE,JDS,JDE,KDS,KDE &
404	& ,IMS,IME,JMS,JME,KMS,KME &
405	& ,ITS,ITE,JTS,JTE,KTS,KTE)
406	!
407	!*** SAVE THE NEW TKE AND MIXING LENGTH.
408	!
409	DO K=KTS,KTE
410	KFLIP=KTE+1-K
411	Q2K(KFLIP)=AMAX1(Q2K(KFLIP),EPSQ2L)
412	TKE(I,K,J)=0.5*Q2K(KFLIP)
413	IF(KFLIP/=KTE)EL_MYJ(I,K,J)=EL(KFLIP) ! EL IS NOT DEFINED AT KTE
414	ENDDO
415	!
416	ENDDO
417	!
418	!----------------------------------------------------------------------
419	ENDDO setup_integration
420	!----------------------------------------------------------------------
421	!
422	!*** CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
423	!
424	DO J=JTS,JTE
425	DO I=ITS,ITE
426	PSFC=PINT(I,LOWLYR(I,J),J)
427	THSK(I,J)=TSK(I,J)(1.E5/PSFC)*CAPA
428	ENDDO
429	ENDDO
430	!
431	!----------------------------------------------------------------------
432	!
433	!----------------------------------------------------------------------
434	main_integration: DO J=JTS,JTE
435	!----------------------------------------------------------------------
436	!
437	DO I=ITS,ITE
438	!
439	!*** FILL 1-D VERTICAL ARRAYS
440	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
441	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
442	!
443	DO K=KTE,KTS,-1
444	KFLIP=KTE+1-K
445	THEK(K)=THE(I,KFLIP,J)
446	RATIOMX=QV(I,KFLIP,J)
447	QK(K)=RATIOMX/(1.+RATIOMX)
448	CWMK(K)=CWM(I,KFLIP,J)
449	ZHK(K)=ZINT(I,K,J)
450	RHOK(K,I)=PMID(I,KFLIP,J)/(R_DT(I,KFLIP,J) &
451	& (1.+P608*QK(K)-CWMK(K)))
452	ENDDO
453	!
454	!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
455	!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. THESE COEFFICIENTS
456	!*** ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
457	!
458	DO K=KTS,KTE-1
459	AKHK(K)=AKH(I,K,J)0.5(RHOK(K,I)+RHOK(K+1,I))
460	ENDDO
461	!
462	ZHK(KTE+1)=ZINT(I,KTE+1,J)
463	!
464	SEAMASK=XLAND(I,J)-1.
465	THZ0(I,J)=(1.-SEAMASK)THSK(I,J)+SEAMASKTHZ0(I,J)
466	!
467	LLOW=LOWLYR(I,J)
468	AKHS_DENS=AKHS(I,J)*RHOK(KTE+1-LLOW,I)
469	!
470	IF(SEAMASK<0.5)THEN
471	QFC1=XLVCHKLOWQ(I,J)AKHS_DENS
472	!
473	IF(SNOW(I,J)>0..OR.SICE(I,J)>0.5)THEN
474	QFC1=QFC1*RLIVWV
475	ENDIF
476	!
477	IF(QFC1>0.)THEN
478	QLOW=QK(KTE+1-LLOW)
479	QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
480	ENDIF
481	!
482	ELSE
483	PSFC=PINT(I,LOWLYR(I,J),J)
484	EXNSFC=(1.E5/PSFC)**CAPA
485	QSFC(I,J)=PQ0SEA/PSFC &
486	& EXP(A2(THSK(I,J)-A3EXNSFC)/(THSK(I,J)-A4EXNSFC))
487	ENDIF
488	!
489	QZ0 (I,J)=(1.-SEAMASK)QSFC(I,J)+SEAMASKQZ0 (I,J)
490	!
491	!*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
492	!
493	LMH=KTE-LOWLYR(I,J)+1
494	!
495	!----------------------------------------------------------------------
496	!*** CARRY OUT THE VERTICAL DIFFUSION OF
497	!*** TEMPERATURE AND WATER VAPOR
498	!----------------------------------------------------------------------
499	!
500	CALL VDIFH(DTDIF,LMH,THZ0(I,J),QZ0(I,J) &
501	& ,AKHS_DENS,CHKLOWQ(I,J),CT(I,J) &
502	& ,THEK,QK,CWMK,AKHK,ZHK,RHOK(KTS,I) &
503	& ,IDS,IDE,JDS,JDE,KDS,KDE &
504	& ,IMS,IME,JMS,JME,KMS,KME &
505	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
506	!----------------------------------------------------------------------
507	!***
508	!*** COMPUTE PRIMARY VARIABLE TENDENCIES
509	!***
510	DO K=KTS,KTE
511	KFLIP=KTE+1-K
512	THOLD=TH(I,K,J)
513	THNEW=THEK(KFLIP)+CWMK(KFLIP)ELOCPAPE(I,K,J)
514	DTDT=(THNEW-THOLD)*RDTTURBL
515	QOLD=QV(I,K,J)/(1.+QV(I,K,J))
516	DQDT=(QK(KFLIP)-QOLD)*RDTTURBL
517	DCDT=(CWMK(KFLIP)-CWM(I,K,J))*RDTTURBL
518	!
519	RTHBLTEN(I,K,J)=DTDT
520	RQVBLTEN(I,K,J)=DQDT/(1.-QK(KFLIP))**2
521	RQCBLTEN(I,K,J)=DCDT
522	ENDDO
523	!
524	!*** Begin debugging
525	! IF(I==IMD.AND.J==JMD)THEN
526	! PRINT_DIAG=0
527	! ELSE
528	! PRINT_DIAG=0
529	! ENDIF
530	! IF(I==227.AND.J==363)PRINT_DIAG=0
531	!*** End debugging
532	!
533	PSFC=.01*PINT(I,LOWLYR(I,J),J)
534	ZSL_DIAG=0.5*DZ(I,1,J)
535	!
536	!*** Begin debugging
537	! IF(PRINT_DIAG==1)THEN
538	!
539	! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
540	! '{turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
541	! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
542	! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
543	! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
544	! '{turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
545	! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
546	! , 2.tke_myj(i,1,j), akhs(i,j), akhs(i,j)ZSL_diag
547	! write(6,"(a)") &
548	! '{turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, S2, EL, Q**2, Akh, EXCH_h, Dz, Dp'
549	! do k=kts,kte/2
550	! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
551	! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
552	! '{turb5 ', k, .01pmid(i,k,j),.01pint(i,k,j), T(i,k,j)-273.15 &
553	! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), S2(KFLIP) &
554	! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
555	! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
556	! enddo
557	!
558	! ELSEIF(PRINT_DIAG==2)THEN
559	!
560	! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
561	! '}turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
562	! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
563	! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
564	! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
565	! '}turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
566	! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
567	! , 2.tke_myj(i,1,j), akhs(i,j), akhs(i,j)ZSL_diag
568	! write(6,"(a)") &
569	! '}turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, S2, EL, Q**2, Akh, EXCH_h, Dz, Dp'
570	! do k=kts,kte/2
571	! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
572	! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
573	! '}turb5 ', k, .01pmid(i,k,j),.01pint(i,k,j), T(i,k,j)-273.15 &
574	! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), S2(KFLIP) &
575	! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
576	! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
577	! enddo
578	! ENDIF
579	!*** End debugging
580	!
581	!----------------------------------------------------------------------
582	ENDDO
583	!----------------------------------------------------------------------
584	DO I=ITS,ITE
585	!
586	!*** FILL 1-D VERTICAL ARRAYS
587	!*** AND FLIP DIRECTION SINCE MYJ SCHEME
588	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
589	!
590	DO K=KTS,KTE-1
591	AKMK(K)=AKM(I,K,J)
592	AKMK(K)=AKMK(K)(RHOK(K,I)+RHOK(K+1,I))0.5
593	ENDDO
594	!
595	LLOW=LOWLYR(I,J)
596	AKMS_DENS=AKMS(I,J)*RHOK(KTE+1-LLOW,I)
597	!
598	DO K=KTE,KTS,-1
599	KFLIP=KTE+1-K
600	UK(K)=U(I,KFLIP,J)
601	VK(K)=V(I,KFLIP,J)
602	ZHK(K)=ZINT(I,K,J)
603	ENDDO
604	ZHK(KTE+1)=ZINT(I,KTE+1,J)
605	!
606	!----------------------------------------------------------------------
607	!*** CARRY OUT THE VERTICAL DIFFUSION OF
608	!*** VELOCITY COMPONENTS
609	!----------------------------------------------------------------------
610	!
611	CALL VDIFV(LMH,DTDIF,UZ0(I,J),VZ0(I,J) &
612	& ,AKMS_DENS,UK,VK,AKMK,ZHK,RHOK(KTS,I) &
613	& ,IDS,IDE,JDS,JDE,KDS,KDE &
614	& ,IMS,IME,JMS,JME,KMS,KME &
615	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
616	!
617	!----------------------------------------------------------------------
618	!***
619	!*** COMPUTE PRIMARY VARIABLE TENDENCIES
620	!***
621	DO K=KTS,KTE
622	KFLIP=KTE+1-K
623	DUDT=(UK(KFLIP)-U(I,K,J))*RDTTURBL
624	DVDT=(VK(KFLIP)-V(I,K,J))*RDTTURBL
625	RUBLTEN(I,K,J)=DUDT
626	RVBLTEN(I,K,J)=DVDT
627	ENDDO
628	!
629	ENDDO
630	!----------------------------------------------------------------------
631	!
632	ENDDO main_integration
633	!
634	!----------------------------------------------------------------------
635	!
636	END SUBROUTINE QNSEPBL
637	!
638	!----------------------------------------------------------------------
639	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
640	!----------------------------------------------------------------------
641	SUBROUTINE MIXLEN &
642	!----------------------------------------------------------------------
643	! ******************************************************************
644	! * *
645	! * LEVEL 2.5 MIXING LENGTH *
646	! * *
647	! ******************************************************************
648	!
649	&(LMH,U,V,T,THE,Q,CWM,Q2,Z,USTAR,CORF &
650	&,S2,GH,RI,EL,PBLH,LPBL,LMXL,CT &
651	&,IDS,IDE,JDS,JDE,KDS,KDE &
652	&,IMS,IME,JMS,JME,KMS,KME &
653	&,ITS,ITE,JTS,JTE,KTS,KTE)
654	!----------------------------------------------------------------------
655	!
656	IMPLICIT NONE
657	!
658	!----------------------------------------------------------------------
659	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
660	& ,IMS,IME,JMS,JME,KMS,KME &
661	& ,ITS,ITE,JTS,JTE,KTS,KTE
662	!
663	INTEGER,INTENT(IN) :: LMH
664	!
665	INTEGER,INTENT(OUT) :: LMXL,LPBL
666	!
667	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: CWM,Q,Q2,T,THE,U,V
668	!
669	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
670	!
671	REAL,INTENT(OUT) :: PBLH
672	!
673	REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: EL,RI,GH,S2
674	!
675	REAL,INTENT(INOUT) :: CT
676	!
677	REAL,INTENT(IN) :: CORF,USTAR
678	!----------------------------------------------------------------------
679	!***
680	!*** LOCAL VARIABLES
681	!***
682	INTEGER :: K,LPBLM
683	!
684	REAL :: A,ADEN,B,BDEN,AUBR,BUBR,BLMX,EL0,ELOQ2X,GHL,S2L &
685	& ,QOL2ST,QOL2UN,QDZL,RDZ,SQ,SREL,SZQ,TEM,THM,VKRMZ,RLAMBDA &
686	& ,RLB,RLN,F
687	!
688	REAL,DIMENSION(KTS:KTE) :: Q1,EN2
689	!
690	REAL,DIMENSION(KTS:KTE-1) :: DTH,ELM,REL
691	!
692	!----------------------------------------------------------------------
693	!**********************************************************************
694	!--------------FIND THE HEIGHT OF THE PBL-------------------------------
695	LPBL=LMH
696	!
697	DO K=LMH-1,1,-1
698	IF(Q2(K)<=EPSQ2L*FH)THEN
699	LPBL=K
700	GO TO 110
701	ENDIF
702	ENDDO
703	!
704	LPBL=1
705	!
706	!--------------THE HEIGHT OF THE PBL------------------------------------
707	!
708	110 PBLH=Z(LPBL)-Z(LMH+1)
709	!
710	!-----------------------------------------------------------------------
711	DO K=KTS,LMH
712	Q1(K)=0.
713	ENDDO
714	!
715	DO K=1,LMH-1
716	DTH(K)=THE(K)-THE(K+1)
717	ENDDO
718	!
719	DO K=LMH-2,1,-1
720	IF(DTH(K)>0..AND.DTH(K+1)<=0.)THEN
721	DTH(K)=DTH(K)+CT
722	EXIT
723	ENDIF
724	ENDDO
725	!
726	CT=0.
727	!----------------------------------------------------------------------
728	!*** COMPUTE LOCAL GRADIENT RICHARDSON NUMBER
729	!----------------------------------------------------------------------
730	DO K=KTS,LMH-1
731	RDZ=2./(Z(K)-Z(K+2))
732	S2L=((U(K)-U(K+1))2+(V(K)-V(K+1))2)RDZRDZ ! S**2
733	S2L=MAX(S2L,EPSGM)
734	S2(K)=S2L
735	!
736	TEM=(T(K)+T(K+1))*0.5
737	THM=(THE(K)+THE(K+1))*0.5
738	!
739	A=THM*P608
740	B=(ELOCP/TEM-1.-P608)*THM
741	!
742	GHL=(DTH(K)((Q(K)+Q(K+1)+CWM(K)+CWM(K+1))(0.5*P608)+1.) &
743	& +(Q(K)-Q(K+1)+CWM(K)-CWM(K+1))*A &
744	& +(CWM(K)-CWM(K+1))B)RDZ ! dTheta/dz
745	!
746	IF(ABS(GHL)<=EPSGH)GHL=EPSGH
747	!
748	EN2(K)=GHLG/THM ! N*2
749	!
750	GH(K)=GHL
751	RI(K)=EN2(K)/S2L
752	ENDDO
753	!
754	!----------------------------------------------------------------------
755	!*** FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
756	!----------------------------------------------------------------------
757	!
758	LMXL=LMH
759	!
760	DO K=KTS,LMH-1
761	S2L=S2(K)
762	GHL=GH(K)
763	!
764	IF(GHL>=EPSGH)THEN
765	IF(S2L/GHL<=REQU)THEN
766	ELM(K)=EPSL
767	LMXL=K
768	ELSE
769	AUBR=(AUBMS2L+AUBHGHL)*GHL
770	BUBR= BUBMS2L+BUBHGHL
771	QOL2ST=(-0.5BUBR+SQRT(BUBRBUBR0.25-AUBRCUBR))*RCUBR
772	ELOQ2X=1./QOL2ST
773	ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
774	ENDIF
775	ELSE
776	ADEN=(ADNMS2L+ADNHGHL)*GHL
777	BDEN= BDNMS2L+BDNHGHL
778	QOL2UN=-0.5BDEN+SQRT(BDENBDEN*0.25-ADEN)
779	ELOQ2X=1./(QOL2UN+EPSRU) ! repsr1/qol2un
780	ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
781	ENDIF
782	ENDDO
783	!
784	IF(ELM(LMH-1)==EPSL)LMXL=LMH
785	!
786	!----------------------------------------------------------------------
787	!*** THE HEIGHT OF THE MIXED LAYER
788	!----------------------------------------------------------------------
789	!
790	BLMX=Z(LMXL)-Z(LMH+1)
791	!
792	!----------------------------------------------------------------------
793	DO K=LPBL,LMH
794	Q1(K)=SQRT(Q2(K))
795	ENDDO
796	!----------------------------------------------------------------------
797	SZQ=0.
798	SQ =0.
799	!
800	DO K=KTS,LMH-1
801	QDZL=(Q1(K)+Q1(K+1))*(Z(K+1)-Z(K+2))
802	SZQ=(Z(K+1)+Z(K+2)-Z(LMH+1)-Z(LMH+1))*QDZL+SZQ
803	SQ=QDZL+SQ
804	ENDDO
805	!
806	!----------------------------------------------------------------------
807	!*** COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
808	!----------------------------------------------------------------------
809	!
810	EL0=MIN(ALPHSZQ0.5/SQ,EL0MAX)
811	EL0=MAX(EL0 ,EL0MIN)
812	!
813	!----------------------------------------------------------------------
814	!*** ABOVE THE PBL TOP
815	!----------------------------------------------------------------------
816	!
817	LPBLM=MAX(LPBL-1,1)
818	!
819	DO K=KTS,LPBLM
820	EL(K)=MIN((Z(K)-Z(K+2))*ELFC,ELM(K))
821	REL(K)=EL(K)/ELM(K)
822	ENDDO
823	!
824	!----------------------------------------------------------------------
825	!*** INSIDE THE PBL
826	!----------------------------------------------------------------------
827	!
828	IF(LPBL<LMH)THEN
829	DO K=LPBL,LMH-1
830	VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
831	EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))
832	REL(K)=EL(K)/ELM(K)
833	ENDDO
834	ENDIF
835	!
836	DO K=LPBL+1,LMH-2
837	SREL=MIN(((REL(K-1)+REL(K+1))0.5+REL(K))0.5,REL(K))
838	EL(K)=MAX(SREL*ELM(K),EPSL)
839	ENDDO
840	!
841	!----------------------------------------------------------------------
842	!*** MIXING LENGTH FOR THE QNSE MODEL IN STABLE CASE
843	!----------------------------------------------------------------------
844	!
845	F=MAX(CORF,EPS1)
846	RLAMBDA=F/(BLCKDR*USTAR)
847	DO K=KTS,LMH-1
848	IF(EN2(K)>=0.0)THEN ! Stable case
849	VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
850	RLB=RLAMBDA+1./VKRMZ
851	RLN=SQRT(2.*EN2(K)/Q2(K))/CN
852	! EL(K)=MIN(1./(RLB+RLN),ELM(K))
853	EL(K)=1./(RLB+RLN)
854	ENDIF
855	ENDDO
856	!
857	!----------------------------------------------------------------------
858	END SUBROUTINE MIXLEN
859	!----------------------------------------------------------------------
860	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
861	!----------------------------------------------------------------------
862	SUBROUTINE PRODQ2 &
863	!----------------------------------------------------------------------
864	! ******************************************************************
865	! * *
866	! * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
867	! * *
868	! ******************************************************************
869	!
870	&(LMH,DTTURBL,USTAR,S2,RI,Q2,EL,Z,AKM,AKH &
871	&,IDS,IDE,JDS,JDE,KDS,KDE &
872	&,IMS,IME,JMS,JME,KMS,KME &
873	&,ITS,ITE,JTS,JTE,KTS,KTE)
874	!----------------------------------------------------------------------
875	!
876	IMPLICIT NONE
877	!
878	!----------------------------------------------------------------------
879	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
880	& ,IMS,IME,JMS,JME,KMS,KME &
881	& ,ITS,ITE,JTS,JTE,KTS,KTE
882	!
883	INTEGER,INTENT(IN) :: LMH
884	!
885	REAL,INTENT(IN) :: DTTURBL,USTAR
886	!
887	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: S2,RI,AKM,AKH,EL
888	!
889	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
890	!
891	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
892	!----------------------------------------------------------------------
893	!***
894	!*** LOCAL VARIABLES
895	!***
896	INTEGER :: K
897	!
898	REAL :: S2L,Q2L,DELTAZ,AKML,AKHL,EN2,PR,BPR,DIS,RC02
899	!
900	!----------------------------------------------------------------------
901	!**********************************************************************
902	!----------------------------------------------------------------------
903	!
904	RC02=2.0/(C0*C0)
905	main_integration: DO K=1,LMH-1
906	S2L=S2(K)
907	Q2L=Q2(K)
908	DELTAZ=0.5*(Z(K)-Z(K+2))
909	AKML=AKM(K)*DELTAZ
910	AKHL=AKH(K)*DELTAZ
911	EN2=RI(K)S2L !N*2
912	!
913	!*** TURBULENCE PRODUCTION TERM
914	!
915	PR=AKML*S2L
916	!
917	!*** BUOYANCY PRODUCTION
918	!
919	BPR=AKHL*EN2
920	!
921	!*** DISSIPATION
922	!
923	DIS=CEPS(0.5Q2L)**1.5/EL(K)
924	!
925	Q2L=Q2L+2.0(PR-BPR-DIS)DTTURBL
926	Q2(K)=AMAX1(Q2L,EPSQ2L)
927	!----------------------------------------------------------------------
928	!*** END OF PRODUCTION/DISSIPATION LOOP
929	!----------------------------------------------------------------------
930	!
931	ENDDO main_integration
932	!
933	!----------------------------------------------------------------------
934	!*** LOWER BOUNDARY CONDITION FOR Q2
935	!----------------------------------------------------------------------
936	!
937	Q2(LMH)=AMAX1(RC02USTARUSTAR,EPSQ2L)
938	!----------------------------------------------------------------------
939	!
940	END SUBROUTINE PRODQ2
941	!
942	!----------------------------------------------------------------------
943	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
944	!----------------------------------------------------------------------
945	SUBROUTINE DIFCOF &
946	! ******************************************************************
947	! * *
948	! * DIFFUSION COEFFICIENTS KM, KH BASED ON THE QNSE THEORY *
949	! * *
950	! ******************************************************************
951	&(LMH,EL,RI,Q2,Z,AKM,AKH &
952	&,IDS,IDE,JDS,JDE,KDS,KDE &
953	&,IMS,IME,JMS,JME,KMS,KME &
954	&,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
955	!----------------------------------------------------------------------
956	!
957	IMPLICIT NONE
958	!
959	!----------------------------------------------------------------------
960	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
961	& ,IMS,IME,JMS,JME,KMS,KME &
962	& ,ITS,ITE,JTS,JTE,KTS,KTE
963	!
964	INTEGER,INTENT(IN) :: LMH
965	!
966	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: Q2
967	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL,RI
968	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
969	!
970	REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: AKH,AKM
971	!----------------------------------------------------------------------
972	!***
973	!*** LOCAL VARIABLES
974	!***
975	INTEGER :: K
976	!
977	REAL :: AK0,ALPHAM,ALPHAH,RIL,RIL2,ARIL,ARIL2,ARIL4,ELL,Q1L,RDZ &
978	& ,AK0DZ,AKMIN
979	!
980	!*** Begin debugging
981	INTEGER,INTENT(IN) :: PRINT_DIAG
982	! REAL :: D2Tmin
983	!*** End debugging
984	!
985	!----------------------------------------------------------------------
986	!**********************************************************************
987	!----------------------------------------------------------------------
988	!
989	DO K=1,LMH-1
990	ELL=EL(K)
991	Q1L=SQRT(0.5*Q2(K)) !Note that Q1L is SQRT(TKE)
992	RIL=RI(K)
993	AK0=C0ELLQ1L !KM in neutral case
994	!
995	!----------------------------------------------------------------------
996	!*** STABILITY FUNCTIONS ALPHAM AND ALPHAH
997	!----------------------------------------------------------------------
998	!
999	!!! UNSTABLE CASE
1000	!
1001	IF(RIL<=0) THEN
1002	ARIL=MIN(ABS(RIL),2.*ARIMIN)
1003	ARIL2=ARIL*ARIL
1004	ARIL4=ARIL2*ARIL2
1005	ALPHAM=1.0+BM1ARIL+BM2ARIL2
1006	ALPHAH=AH1+BH1ARIL+BH2ARIL2+BH3*ARIL4
1007	!
1008	!!! STABLE CASE
1009	!
1010	ELSE
1011	RIL2=RIL*RIL
1012	ALPHAM=(1.0+AM1RIL2)/(1.0+AM2RIL+AM3*RIL2)
1013	ALPHAH=(AH1+AH2RIL+AH3RIL2)/(1.0+AH4RIL+AH5RIL2)
1014	ENDIF
1015	!
1016	!-----------------------------------------------------------------------
1017	!*** END OF STABILITY FUNCTIONS COMPUTATIONS
1018	!-----------------------------------------------------------------------
1019	!
1020	!----------------------------------------------------------------------
1021	!*** DIFFUSION COEFFICIENTS
1022	!----------------------------------------------------------------------
1023	!
1024	RDZ=2./(Z(K)-Z(K+2))
1025	AK0DZ=AK0*RDZ
1026	AKMIN=EPSKM*RDZ
1027	AKM(K)=MAX(AK0DZ*ALPHAM,AKMIN)
1028	AKH(K)=MAX(AK0DZ*ALPHAH,AKMIN)
1029	!----------------------------------------------------------------------
1030	ENDDO
1031	!----------------------------------------------------------------------
1032	!
1033	END SUBROUTINE DIFCOF
1034	!
1035	!----------------------------------------------------------------------
1036	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1037	!----------------------------------------------------------------------
1038	SUBROUTINE VDIFQ &
1039	! ******************************************************************
1040	! * *
1041	! * VERTICAL DIFFUSION OF Q2 (TKE) *
1042	! * *
1043	! ******************************************************************
1044	&(LMH,DTDIF,Q2,EL,Z &
1045	&,IDS,IDE,JDS,JDE,KDS,KDE &
1046	&,IMS,IME,JMS,JME,KMS,KME &
1047	&,ITS,ITE,JTS,JTE,KTS,KTE)
1048	!----------------------------------------------------------------------
1049	!
1050	IMPLICIT NONE
1051	!
1052	!----------------------------------------------------------------------
1053	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1054	& ,IMS,IME,JMS,JME,KMS,KME &
1055	& ,ITS,ITE,JTS,JTE,KTS,KTE
1056	!
1057	INTEGER,INTENT(IN) :: LMH
1058	!
1059	REAL,INTENT(IN) :: DTDIF
1060	!
1061	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL
1062	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1063	!
1064	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
1065	!----------------------------------------------------------------------
1066	!***
1067	!*** LOCAL VARIABLES
1068	!***
1069	INTEGER :: K
1070	!
1071	REAL :: ADEN,AKQS,BDEN,BESH,BESM,CDEN,CF,DTOZS,ELL,ELOQ2,ELOQ4 &
1072	& ,ELQDZ,ESH,ESM,ESQHF,GHL,GML,Q1L,RDEN,RDZ
1073	!
1074	REAL,DIMENSION(KTS:KTE-2) :: AKQ,CM,CR,DTOZ,RSQ2
1075	!----------------------------------------------------------------------
1076	!**********************************************************************
1077	!----------------------------------------------------------------------
1078	!***
1079	!*** VERTICAL TURBULENT DIFFUSION
1080	!***
1081	!----------------------------------------------------------------------
1082	ESQHF=0.5*ESQ
1083	!
1084	DO K=KTS,LMH-2
1085	DTOZ(K)=(DTDIF+DTDIF)/(Z(K)-Z(K+2))
1086	AKQ(K)=SQRT((Q2(K)+Q2(K+1))0.5)(EL(K)+EL(K+1))*ESQHF &
1087	& /(Z(K+1)-Z(K+2))
1088	CR(K)=-DTOZ(K)*AKQ(K)
1089	ENDDO
1090	!
1091	CM(1)=DTOZ(1)*AKQ(1)+1.
1092	RSQ2(1)=Q2(1)
1093	!
1094	DO K=KTS+1,LMH-2
1095	CF=-DTOZ(K)*AKQ(K-1)/CM(K-1)
1096	CM(K)=-CR(K-1)CF+(AKQ(K-1)+AKQ(K))DTOZ(K)+1.
1097	RSQ2(K)=-RSQ2(K-1)*CF+Q2(K)
1098	ENDDO
1099	!
1100	DTOZS=(DTDIF+DTDIF)/(Z(LMH-1)-Z(LMH+1))
1101	AKQS=SQRT((Q2(LMH-1)+Q2(LMH))0.5)(EL(LMH-1)+ELZ0)*ESQHF &
1102	& /(Z(LMH)-Z(LMH+1))
1103	!
1104	CF=-DTOZS*AKQ(LMH-2)/CM(LMH-2)
1105	!
1106	Q2(LMH-1)=(DTOZSAKQSQ2(LMH)-RSQ2(LMH-2)*CF+Q2(LMH-1)) &
1107	& /((AKQ(LMH-2)+AKQS)DTOZS-CR(LMH-2)CF+1.)
1108	!
1109	DO K=LMH-2,KTS,-1
1110	Q2(K)=(-CR(K)*Q2(K+1)+RSQ2(K))/CM(K)
1111	ENDDO
1112	!----------------------------------------------------------------------
1113	!
1114	END SUBROUTINE VDIFQ
1115	!
1116	!----------------------------------------------------------------------
1117	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1118	!---------------------------------------------------------------------
1119	SUBROUTINE VDIFH(DTDIF,LMH,THZ0,QZ0,RKHS,CHKLOWQ,CT &
1120	& ,THE,Q,CWM,RKH,Z,RHO &
1121	& ,IDS,IDE,JDS,JDE,KDS,KDE &
1122	& ,IMS,IME,JMS,JME,KMS,KME &
1123	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1124	! ***************************************************************
1125	! * *
1126	! * VERTICAL DIFFUSION OF MASS VARIABLES *
1127	! * *
1128	! ***************************************************************
1129	!---------------------------------------------------------------------
1130	!
1131	IMPLICIT NONE
1132	!
1133	!---------------------------------------------------------------------
1134	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1135	& ,IMS,IME,JMS,JME,KMS,KME &
1136	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1137	!
1138	INTEGER,INTENT(IN) :: LMH
1139	!
1140	REAL,INTENT(IN) :: CHKLOWQ,CT,DTDIF,QZ0,RKHS,THZ0
1141	!
1142	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKH
1143	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1144	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1145	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: CWM,Q,THE
1146	!
1147	!----------------------------------------------------------------------
1148	!***
1149	!*** LOCAL VARIABLES
1150	!***
1151	INTEGER :: K
1152	!
1153	REAL :: CF,CMB,CMCB,CMQB,CMTB,CTHF,DTOZL,DTOZS &
1154	& ,RCML,RKHH,RKQS,RSCB,RSQB,RSTB
1155	!
1156	REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RKCT,RSC,RSQ,RST
1157	!
1158	!----------------------------------------------------------------------
1159	!**********************************************************************
1160	!----------------------------------------------------------------------
1161	CTHF=0.5*CT
1162	!
1163	DO K=KTS,LMH-1
1164	DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1165	CR(K)=-DTOZ(K)*RKH(K)
1166	RKCT(K)=RKH(K)(Z(K)-Z(K+2))CTHF
1167	ENDDO
1168	!
1169	CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHO(KTS)
1170	!----------------------------------------------------------------------
1171	RST(KTS)=-RKCT(KTS)*DTOZ(KTS) &
1172	& +THE(KTS)*RHO(KTS)
1173	RSQ(KTS)=Q(KTS) *RHO(KTS)
1174	RSC(KTS)=CWM(KTS)*RHO(KTS)
1175	!----------------------------------------------------------------------
1176	DO K=KTS+1,LMH-1
1177	DTOZL=DTOZ(K)
1178	CF=-DTOZL*RKH(K-1)/CM(K-1)
1179	CM(K)=-CR(K-1)CF+(RKH(K-1)+RKH(K))DTOZL+RHO(K)
1180	RST(K)=-RST(K-1)CF+(RKCT(K-1)-RKCT(K))DTOZL+THE(K)*RHO(K)
1181	RSQ(K)=-RSQ(K-1)CF+Q(K) RHO(K)
1182	RSC(K)=-RSC(K-1)CF+CWM(K)RHO(K)
1183	ENDDO
1184	!
1185	DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1186	RKHH=RKH(LMH-1)
1187	!
1188	CF=-DTOZS*RKHH/CM(LMH-1)
1189	RKQS=RKHS*CHKLOWQ
1190	!
1191	CMB=CR(LMH-1)*CF
1192	CMTB=-CMB+(RKHH+RKHS)*DTOZS+RHO(LMH)
1193	CMQB=-CMB+(RKHH+RKQS)*DTOZS+RHO(LMH)
1194	CMCB=-CMB+(RKHH )*DTOZS+RHO(LMH)
1195	!
1196	RSTB=-RST(LMH-1)CF+RKCT(LMH-1)DTOZS+THE(LMH)*RHO(LMH)
1197	RSQB=-RSQ(LMH-1)CF+Q(LMH) RHO(LMH)
1198	RSCB=-RSC(LMH-1)CF+CWM(LMH)RHO(LMH)
1199	!----------------------------------------------------------------------
1200	THE(LMH)=(DTOZSRKHSTHZ0+RSTB)/CMTB
1201	Q(LMH) =(DTOZSRKQSQZ0 +RSQB)/CMQB
1202	CWM(LMH)=( RSCB)/CMCB
1203	!----------------------------------------------------------------------
1204	DO K=LMH-1,KTS,-1
1205	RCML=1./CM(K)
1206	THE(K)=(-CR(K)THE(K+1)+RST(K))RCML
1207	Q(K) =(-CR(K)* Q(K+1)+RSQ(K))*RCML
1208	CWM(K)=(-CR(K)CWM(K+1)+RSC(K))RCML
1209	ENDDO
1210	!----------------------------------------------------------------------
1211	!
1212	END SUBROUTINE VDIFH
1213	!
1214	!---------------------------------------------------------------------
1215	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1216	!---------------------------------------------------------------------
1217	SUBROUTINE VDIFV(LMH,DTDIF,UZ0,VZ0,RKMS,U,V,RKM,Z,RHO &
1218	& ,IDS,IDE,JDS,JDE,KDS,KDE &
1219	& ,IMS,IME,JMS,JME,KMS,KME &
1220	,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1221	! ***************************************************************
1222	! * *
1223	! * VERTICAL DIFFUSION OF VELOCITY COMPONENTS *
1224	! * *
1225	! ***************************************************************
1226	!---------------------------------------------------------------------
1227	!
1228	IMPLICIT NONE
1229	!
1230	!---------------------------------------------------------------------
1231	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1232	& ,IMS,IME,JMS,JME,KMS,KME &
1233	& ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1234	!
1235	INTEGER,INTENT(IN) :: LMH
1236	!
1237	REAL,INTENT(IN) :: RKMS,DTDIF,UZ0,VZ0
1238	!
1239	REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKM
1240	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1241	REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1242	!
1243	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: U,V
1244	!----------------------------------------------------------------------
1245	!***
1246	!*** LOCAL VARIABLES
1247	!***
1248	INTEGER :: K
1249	!
1250	REAL :: CF,DTOZAK,DTOZL,DTOZS,RCML,RCMVB,RHOK,RKMH
1251	!
1252	REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RSU,RSV
1253	!----------------------------------------------------------------------
1254	!**********************************************************************
1255	!----------------------------------------------------------------------
1256	DO K=1,LMH-1
1257	DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1258	CR(K)=-DTOZ(K)*RKM(K)
1259	ENDDO
1260	!
1261	RHOK=RHO(1)
1262	CM(1)=DTOZ(1)*RKM(1)+RHOK
1263	RSU(1)=U(1)*RHOK
1264	RSV(1)=V(1)*RHOK
1265	!----------------------------------------------------------------------
1266	DO K=2,LMH-1
1267	DTOZL=DTOZ(K)
1268	CF=-DTOZL*RKM(K-1)/CM(K-1)
1269	RHOK=RHO(K)
1270	CM(K)=-CR(K-1)CF+(RKM(K-1)+RKM(K))DTOZL+RHOK
1271	RSU(K)=-RSU(K-1)CF+U(K)RHOK
1272	RSV(K)=-RSV(K-1)CF+V(K)RHOK
1273	ENDDO
1274	!----------------------------------------------------------------------
1275	DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1276	RKMH=RKM(LMH-1)
1277	!
1278	CF=-DTOZS*RKMH/CM(LMH-1)
1279	RHOK=RHO(LMH)
1280	RCMVB=1./((RKMH+RKMS)DTOZS-CR(LMH-1)CF+RHOK)
1281	DTOZAK=DTOZS*RKMS
1282	!----------------------------------------------------------------------
1283	U(LMH)=(DTOZAKUZ0-RSU(LMH-1)CF+U(LMH)RHOK)RCMVB
1284	V(LMH)=(DTOZAKVZ0-RSV(LMH-1)CF+V(LMH)RHOK)RCMVB
1285	!----------------------------------------------------------------------
1286	DO K=LMH-1,1,-1
1287	RCML=1./CM(K)
1288	U(K)=(-CR(K)U(K+1)+RSU(K))RCML
1289	V(K)=(-CR(K)V(K+1)+RSV(K))RCML
1290	ENDDO
1291	!----------------------------------------------------------------------
1292	!
1293	END SUBROUTINE VDIFV
1294	!
1295	!-----------------------------------------------------------------------
1296	!
1297	!=======================================================================
1298	SUBROUTINE QNSEPBLINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1299	& TKE,EXCH_H,RESTART,ALLOWED_TO_READ, &
1300	& IDS,IDE,JDS,JDE,KDS,KDE, &
1301	& IMS,IME,JMS,JME,KMS,KME, &
1302	& ITS,ITE,JTS,JTE,KTS,KTE )
1303	!-----------------------------------------------------------------------
1304	IMPLICIT NONE
1305	!-----------------------------------------------------------------------
1306	LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1307	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1308	& IMS,IME,JMS,JME,KMS,KME, &
1309	& ITS,ITE,JTS,JTE,KTS,KTE
1310
1311	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
1312	& RUBLTEN, &
1313	& RVBLTEN, &
1314	& RTHBLTEN, &
1315	& RQVBLTEN, &
1316	& TKE
1317	INTEGER :: I,J,K,ITF,JTF,KTF
1318	!-----------------------------------------------------------------------
1319	!-----------------------------------------------------------------------
1320
1321	JTF=MIN0(JTE,JDE-1)
1322	KTF=MIN0(KTE,KDE-1)
1323	ITF=MIN0(ITE,IDE-1)
1324
1325	IF(.NOT.RESTART)THEN
1326	DO J=JTS,JTF
1327	DO K=KTS,KTF
1328	DO I=ITS,ITF
1329	TKE(I,K,J)=EPSQ2L
1330	RUBLTEN(I,K,J)=0.
1331	RVBLTEN(I,K,J)=0.
1332	RTHBLTEN(I,K,J)=0.
1333	RQVBLTEN(I,K,J)=0.
1334	EXCH_H(I,K,J)=0.
1335	ENDDO
1336	ENDDO
1337	ENDDO
1338	ENDIF
1339
1340	END SUBROUTINE QNSEPBLINIT
1341	!-----------------------------------------------------------------------
1342	!
1343	END MODULE MODULE_BL_QNSEPBL
1344	!
1345	!-----------------------------------------------------------------------

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

Download in other formats:

Original Format