Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_gfs.F @ 1379

Last change on this file since 1379 was 1, checked in by lfita, 11 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: 41.5 KB

Line
1	!LWRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE module_bl_gfs
4
5	CONTAINS
6
7	!-------------------------------------------------------------------
8	SUBROUTINE BL_GFS(U3D,V3D,TH3D,T3D,QV3D,QC3D,QI3D,P3D,PI3D, &
9	RUBLTEN,RVBLTEN,RTHBLTEN, &
10	RQVBLTEN,RQCBLTEN,RQIBLTEN, &
11	CP,G,ROVCP,R,ROVG,P_QI,P_FIRST_SCALAR, &
12	dz8w,z,PSFC, &
13	UST,PBL,PSIM,PSIH, &
14	HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
15	DT,KPBL2D,EP1,KARMAN, &
16	#if (NMM_CORE==1)
17	DISHEAT, &
18	#endif
19	ids,ide, jds,jde, kds,kde, &
20	ims,ime, jms,jme, kms,kme, &
21	its,ite, jts,jte, kts,kte )
22	!--------------------------------------------------------------------
23	USE MODULE_GFS_MACHINE , ONLY : kind_phys
24	!-------------------------------------------------------------------
25	IMPLICIT NONE
26	!-------------------------------------------------------------------
27	!-- U3D 3D u-velocity interpolated to theta points (m/s)
28	!-- V3D 3D v-velocity interpolated to theta points (m/s)
29	!-- TH3D 3D potential temperature (K)
30	!-- T3D temperature (K)
31	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
32	!-- QC3D 3D cloud mixing ratio (Kg/Kg)
33	!-- QI3D 3D ice mixing ratio (Kg/Kg)
34	!-- P3D 3D pressure (Pa)
35	!-- PI3D 3D exner function (dimensionless)
36	!-- rr3D 3D dry air density (kg/m^3)
37	!-- RUBLTEN U tendency due to
38	! PBL parameterization (m/s^2)
39	!-- RVBLTEN V tendency due to
40	! PBL parameterization (m/s^2)
41	!-- RTHBLTEN Theta tendency due to
42	! PBL parameterization (K/s)
43	!-- RQVBLTEN Qv tendency due to
44	! PBL parameterization (kg/kg/s)
45	!-- RQCBLTEN Qc tendency due to
46	! PBL parameterization (kg/kg/s)
47	!-- RQIBLTEN Qi tendency due to
48	! PBL parameterization (kg/kg/s)
49	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
50	!-- G acceleration due to gravity (m/s^2)
51	!-- ROVCP R/CP
52	!-- R gas constant for dry air (J/kg/K)
53	!-- ROVG R/G
54	!-- P_QI species index for cloud ice
55	!-- dz8w dz between full levels (m)
56	!-- z height above sea level (m)
57	!-- PSFC pressure at the surface (Pa)
58	!-- UST u* in similarity theory (m/s)
59	!-- PBL PBL height (m)
60	!-- PSIM similarity stability function for momentum
61	!-- PSIH similarity stability function for heat
62	!-- HFX upward heat flux at the surface (W/m^2)
63	!-- QFX upward moisture flux at the surface (kg/m^2/s)
64	!-- TSK surface temperature (K)
65	!-- GZ1OZ0 log(z/z0) where z0 is roughness length
66	!-- WSPD wind speed at lowest model level (m/s)
67	!-- BR bulk Richardson number in surface layer
68	!-- DT time step (s)
69	!-- rvovrd R_v divided by R_d (dimensionless)
70	!-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
71	!-- KARMAN Von Karman constant
72	!-- ids start index for i in domain
73	!-- ide end index for i in domain
74	!-- jds start index for j in domain
75	!-- jde end index for j in domain
76	!-- kds start index for k in domain
77	!-- kde end index for k in domain
78	!-- ims start index for i in memory
79	!-- ime end index for i in memory
80	!-- jms start index for j in memory
81	!-- jme end index for j in memory
82	!-- kms start index for k in memory
83	!-- kme end index for k in memory
84	!-- its start index for i in tile
85	!-- ite end index for i in tile
86	!-- jts start index for j in tile
87	!-- jte end index for j in tile
88	!-- kts start index for k in tile
89	!-- kte end index for k in tile
90	!-------------------------------------------------------------------
91
92	#if (NMM_CORE==1)
93	LOGICAL , INTENT(IN):: DISHEAT ! gopal's doing
94	#endif
95
96	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
97	ims,ime, jms,jme, kms,kme, &
98	its,ite, jts,jte, kts,kte, &
99	P_QI,P_FIRST_SCALAR
100
101	REAL, INTENT(IN) :: &
102	CP, &
103	DT, &
104	EP1, &
105	G, &
106	KARMAN, &
107	R, &
108	ROVCP, &
109	ROVG
110
111	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
112	DZ8W, &
113	P3D, &
114	PI3D, &
115	QC3D, &
116	QI3D, &
117	QV3D, &
118	T3D, &
119	TH3D, &
120	U3D, &
121	V3D, &
122	Z
123
124
125	REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
126	RTHBLTEN, &
127	RQCBLTEN, &
128	RQIBLTEN, &
129	RQVBLTEN, &
130	RUBLTEN, &
131	RVBLTEN
132
133	REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
134	BR, &
135	GZ1OZ0, &
136	HFX, &
137	PSFC, &
138	PSIM, &
139	PSIH, &
140	QFX, &
141	TSK
142
143	REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
144	PBL, &
145	UST, &
146	WSPD
147
148	INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
149	KPBL2D
150
151
152	!--------------------------- LOCAL VARS ------------------------------
153
154
155	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
156	DEL, &
157	DU, &
158	DV, &
159	PHIL, &
160	PRSL, &
161	PRSLK, &
162	T1, &
163	TAU, &
164	dishx, &
165	U1, &
166	V1
167
168	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
169	PHII, &
170	PRSI
171
172	REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 3) :: &
173	Q1, &
174	RTG
175
176	REAL (kind=kind_phys), DIMENSION(its:ite) :: &
177	DQSFC, &
178	DTSFC, &
179	DUSFC, &
180	DVSFC, &
181	EVAP, &
182	FH, &
183	FM, &
184	HEAT, &
185	HGAMQ, &
186	HGAMT, &
187	HPBL, &
188	PSK, &
189	QSS, &
190	RBSOIL, &
191	RCL, &
192	SPD1, &
193	STRESS, &
194	TSEA
195
196	REAL (kind=kind_phys) :: &
197	CPM, &
198	cpmikj, &
199	DELTIM, &
200	FMTMP, &
201	RRHOX
202
203	INTEGER, DIMENSION( its:ite ) :: &
204	KPBL
205
206	INTEGER :: &
207	I, &
208	IM, &
209	J, &
210	K, &
211	KM, &
212	KTEM, &
213	KTEP, &
214	KX, &
215	L, &
216	NTRAC
217
218	IM=ITE-ITS+1
219	KX=KTE-KTS+1
220	KTEM=KTE-1
221	KTEP=KTE+1
222	NTRAC=2
223	DELTIM=DT
224	IF (P_QI.ge.P_FIRST_SCALAR) NTRAC=3
225
226
227	DO J=jts,jte
228
229	DO i=its,ite
230	RRHOX=(RT3D(I,KTS,J)(1.+EP1*QV3D(I,KTS,J)))/PSFC(I,J)
231	CPM=CP(1.+0.8QV3D(i,kts,j))
232	FMTMP=GZ1OZ0(i,j)-PSIM(i,j)
233	PSK(i)=(PSFC(i,j).00001)*ROVCP
234	FM(i)=FMTMP
235	FH(i)=GZ1OZ0(i,j)-PSIH(i,j)
236	TSEA(i)=TSK(i,j)
237	QSS(i)=QV3D(i,kts,j) ! not used in moninp so set to qv3d for now
238	HEAT(i)=HFX(i,j)/CPM*RRHOX
239	EVAP(i)=QFX(i,j)*RRHOX
240	STRESS(i)=KARMANKARMANWSPD(i,j)WSPD(i,j)/(FMTMPFMTMP)
241	SPD1(i)=WSPD(i,j)
242	PRSI(i,kts)=PSFC(i,j)*.001
243	PHII(I,kts)=0.
244	RCL(i)=1.
245	RBSOIL(I)=BR(i,j)
246	ENDDO
247
248	DO k=kts,kte
249	DO i=its,ite
250	DV(I,K) = 0.
251	DU(I,K) = 0.
252	TAU(I,K) = 0.
253	U1(I,K) = U3D(i,k,j)
254	V1(I,K) = V3D(i,k,j)
255	T1(I,K) = T3D(i,k,j)
256	Q1(I,K,1) = QV3D(i,k,j)/(1.+QV3D(i,k,j))
257	Q1(I,K,2) = QC3D(i,k,j)/(1.+QC3D(i,k,j))
258	PRSL(I,K)=P3D(i,k,j)*.001
259	ENDDO
260	ENDDO
261
262	DO k=kts,kte
263	DO i=its,ite
264	PRSLK(I,K)=(PRSL(i,k).01)*ROVCP
265	ENDDO
266	ENDDO
267
268
269	DO k=kts+1,kte
270	km=k-1
271	DO i=its,ite
272	DEL(i,km)=PRSL(i,km)/ROVG*dz8w(i,km,j)/T3D(i,km,j)
273	PRSI(i,k)=PRSI(i,km)-DEL(i,km)
274	PHII(I,K)=(Z(i,k,j)-Z(i,kts,j))*G
275	PHIL(I,KM)=0.5(Z(i,k,j)+Z(i,km,j)-2.Z(i,kts,j))*G
276	ENDDO
277	ENDDO
278
279	DO i=its,ite
280	DEL(i,kte)=DEL(i,ktem)
281	PRSI(i,ktep)=PRSI(i,kte)-DEL(i,ktem)
282	PHII(I,KTEP)=PHII(I,KTE)+dz8w(i,kte,j)*G
283	PHIL(I,KTE)=PHII(I,KTE)-PHIL(I,KTEM)+PHII(I,KTE)
284	ENDDO
285
286	IF (P_QI.ge.P_FIRST_SCALAR) THEN
287	DO k=kts,kte
288	DO i=its,ite
289	Q1(I,K,3) = QI3D(i,k,j)/(1.+QI3D(i,k,j))
290	ENDDO
291	ENDDO
292	ENDIF
293
294	DO l=1,ntrac
295	DO k=kts,kte
296	DO i=its,ite
297	RTG(I,K,L) = 0.
298	ENDDO
299	ENDDO
300	ENDDO
301
302	CALL MONINP(IM,IM,KX,NTRAC,DV,DU,TAU,RTG,U1,V1,T1,Q1, &
303	PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS, &
304	SPD1,KPBL,PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL, &
305	DELTIM,DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ)
306
307	!============================================================================
308	! ADD IN DISSIPATIVE HEATING .... v*dv. This is Bob's doing
309	!============================================================================
310
311	#if (NMM_CORE==1)
312
313	IF(DISHEAT)THEN
314	DO k=kts,kte
315	DO i=its,ite
316	dishx(i,k)=u1(i,k)du(i,k) + v1(i,k)dv(i,k)
317	cpmikj=CP(1.+0.8QV3D(i,k,j))
318	dishx(i,k)=-dishx(i,k)/cpmikj
319	! IF(k==1)WRITE(0,*)'ADDITIONAL DISSIPATIVE HEATING',tau(i,k),dishx(i,k)
320	tau(i,k)=tau(i,k)+dishx(i,k)
321	ENDDO
322	ENDDO
323	ENDIF
324	#endif
325
326	!=============================================================================
327
328
329	DO k=kts,kte
330	DO i=its,ite
331	RVBLTEN(I,K,J)=DV(I,K)
332	RUBLTEN(I,K,J)=DU(I,K)
333	RTHBLTEN(I,K,J)=TAU(I,K)/PI3D(I,K,J)
334	RQVBLTEN(I,K,J)=RTG(I,K,1)/(1.-Q1(I,K,1))**2
335	RQCBLTEN(I,K,J)=RTG(I,K,2)/(1.-Q1(I,K,2))**2
336	ENDDO
337	ENDDO
338
339	IF (P_QI.ge.P_FIRST_SCALAR) THEN
340	DO k=kts,kte
341	DO i=its,ite
342	RQIBLTEN(I,K,J)=RTG(I,K,3)/(1.-Q1(I,K,3))**2
343	ENDDO
344	ENDDO
345	ENDIF
346
347	DO i=its,ite
348	UST(i,j)=SQRT(STRESS(i))
349	WSPD(i,j)=SQRT(U3D(I,KTS,J)*U3D(I,KTS,J)+ &
350	V3D(I,KTS,J)*V3D(I,KTS,J))+1.E-9
351	PBL(i,j)=HPBL(i)
352	KPBL2D(i,j)=kpbl(i)
353	ENDDO
354
355	ENDDO
356
357
358	END SUBROUTINE BL_GFS
359
360	!===================================================================
361	SUBROUTINE gfsinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
362	RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
363	restart, &
364	allowed_to_read, &
365	ids, ide, jds, jde, kds, kde, &
366	ims, ime, jms, jme, kms, kme, &
367	its, ite, jts, jte, kts, kte )
368	!-------------------------------------------------------------------
369	IMPLICIT NONE
370	!-------------------------------------------------------------------
371	LOGICAL , INTENT(IN) :: allowed_to_read,restart
372	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
373	ims, ime, jms, jme, kms, kme, &
374	its, ite, jts, jte, kts, kte
375	INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
376
377	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
378	RUBLTEN, &
379	RVBLTEN, &
380	RTHBLTEN, &
381	RQVBLTEN, &
382	RQCBLTEN, &
383	RQIBLTEN
384	INTEGER :: i, j, k, itf, jtf, ktf
385
386	jtf=min0(jte,jde-1)
387	ktf=min0(kte,kde-1)
388	itf=min0(ite,ide-1)
389
390	IF(.not.restart)THEN
391	DO j=jts,jtf
392	DO k=kts,ktf
393	DO i=its,itf
394	RUBLTEN(i,k,j)=0.
395	RVBLTEN(i,k,j)=0.
396	RTHBLTEN(i,k,j)=0.
397	RQVBLTEN(i,k,j)=0.
398	RQCBLTEN(i,k,j)=0.
399	ENDDO
400	ENDDO
401	ENDDO
402	ENDIF
403
404	IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
405	DO j=jts,jtf
406	DO k=kts,ktf
407	DO i=its,itf
408	RQIBLTEN(i,k,j)=0.
409	ENDDO
410	ENDDO
411	ENDDO
412	ENDIF
413
414	IF (P_QI .ge. P_FIRST_SCALAR) THEN
415	DO j=jts,jtf
416	DO k=kts,ktf
417	DO i=its,itf
418	RQIBLTEN(i,k,j)=0.
419	ENDDO
420	ENDDO
421	ENDDO
422	ENDIF
423
424
425	END SUBROUTINE gfsinit
426
427	! --------------------------------------------------------------
428	!FPP$ NOCONCUR R
429	SUBROUTINE MONINP(IX,IM,KM,ntrac,DV,DU,TAU,RTG, &
430	& U1,V1,T1,Q1, &
431	& PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS,SPD1,KPBL, &
432	! & PSK,RBSOIL,CD,CH,FM,FH,TSEA,QSS,DPHI,SPD1,KPBL, &
433	& PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,DELTIM, &
434	& DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ)
435	!
436	USE MODULE_GFS_MACHINE, ONLY : kind_phys
437	USE MODULE_GFS_PHYSCONS, grav => con_g, RD => con_RD, CP => con_CP &
438	&, HVAP => con_HVAP, ROG => con_ROG, FV => con_FVirt
439
440	implicit none
441	!
442	! include 'constant.h'
443	!
444	!
445	! Arguments
446	!
447	integer IX, IM, KM, ntrac, KPBL(IM)
448	!
449	real(kind=kind_phys) DELTIM
450	real(kind=kind_phys) DV(IM,KM), DU(IM,KM), &
451	& TAU(IM,KM), RTG(IM,KM,ntrac), &
452	& U1(IX,KM), V1(IX,KM), &
453	& T1(IX,KM), Q1(IX,KM,ntrac), &
454	& PSK(IM), RBSOIL(IM), &
455	! & CD(IM), CH(IM), &
456	& FM(IM), FH(IM), &
457	& TSEA(IM), QSS(IM), &
458	& SPD1(IM), &
459	! & DPHI(IM), SPD1(IM), &
460	& PRSI(IX,KM+1), DEL(IX,KM), &
461	& PRSL(IX,KM), PRSLK(IX,KM), &
462	& PHII(IX,KM+1), PHIL(IX,KM), &
463	& RCL(IM), DUSFC(IM), &
464	& dvsfc(IM), dtsfc(IM), &
465	& DQSFC(IM), HPBL(IM), &
466	& HGAMT(IM), hgamq(IM)
467	!
468	! Locals
469	!
470	integer i,iprt,is,iun,k,kk,kmpbl,lond
471	! real(kind=kind_phys) betaq(IM), betat(IM), betaw(IM), &
472	real(kind=kind_phys) evap(IM), heat(IM), phih(IM), &
473	& phim(IM), rbdn(IM), rbup(IM), &
474	& the1(IM), stress(im), beta(im), &
475	& the1v(IM), thekv(IM), thermal(IM), &
476	& thesv(IM), ustar(IM), wscale(IM)
477	! & thesv(IM), ustar(IM), wscale(IM), zl1(IM)
478	!
479	real(kind=kind_phys) RDZT(IM,KM-1), &
480	& ZI(IM,KM+1), ZL(IM,KM), &
481	& DKU(IM,KM-1), DKT(IM,KM-1), DKO(IM,KM-1), &
482	& AL(IM,KM-1), AD(IM,KM), &
483	& AU(IM,KM-1), A1(IM,KM), &
484	& A2(IM,KM), THETA(IM,KM), &
485	& AT(IM,KM*(ntrac-1))
486	logical pblflg(IM), sfcflg(IM), stable(IM)
487	!
488	real(kind=kind_phys) aphi16, aphi5, bet1, bvf2, &
489	& cfac, conq, cont, conw, &
490	& conwrc, dk, dkmax, dkmin, &
491	& dq1, dsdz2, dsdzq, dsdzt, &
492	& dsig, dt, dthe1, dtodsd, &
493	& dtodsu, dw2, dw2min, g, &
494	& gamcrq, gamcrt, gocp, gor, gravi, &
495	& hol, pfac, prmax, prmin, prinv, &
496	& prnum, qmin, qtend, rbcr, &
497	& rbint, rdt, rdz, rdzt1, &
498	& ri, rimin, rl2, rlam, &
499	& rone, rzero, sfcfrac, &
500	& sflux, shr2, spdk2, sri, &
501	& tem, ti, ttend, tvd, &
502	& tvu, utend, vk, vk2, &
503	& vpert, vtend, xkzo, zfac, &
504	& zfmin, zk, tem1
505	!
506	PARAMETER(g=grav)
507	PARAMETER(GOR=G/RD,GOCP=G/CP)
508	PARAMETER(CONT=1000.CP/G,CONQ=1000.HVAP/G,CONW=1000./G)
509	PARAMETER(RLAM=150.,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.)
510	PARAMETER(DW2MIN=0.0001,DKMIN=1.0,DKMAX=1000.,RIMIN=-100.)
511	PARAMETER(RBCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
512	PARAMETER(QMIN=1.E-8,XKZO=1.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
513	! PARAMETER(GAMCRT=3.,GAMCRQ=2.E-3)
514	PARAMETER(GAMCRT=3.,GAMCRQ=0.)
515	PARAMETER(RZERO=0.,RONE=1.)
516	PARAMETER(IUN=84)
517	!
518	!
519	!-----------------------------------------------------------------------
520	!
521	601 FORMAT(1X,' MONINP LAT LON STEP HOUR ',3I6,F6.1)
522	602 FORMAT(1X,' K',' Z',' T',' TH', &
523	& ' TVH',' Q',' U',' V', &
524	& ' SP')
525	603 FORMAT(1X,I5,8F9.1)
526	604 FORMAT(1X,' SFC',9X,F9.1,18X,F9.1)
527	605 FORMAT(1X,' K ZL SPD2 THEKV THE1V' &
528	& ,' THERMAL RBUP')
529	606 FORMAT(1X,I5,6F8.2)
530	607 FORMAT(1X,' KPBL HPBL FM FH HGAMT', &
531	& ' HGAMQ WS USTAR CD CH')
532	608 FORMAT(1X,I5,9F8.2)
533	609 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2)
534	610 FORMAT(1X,' K PR DKT DKU ',I5,3F8.2,' L2 RI T2', &
535	& ' SR2 ',2F8.2,2E10.2)
536	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
537	! COMPUTE PRELIMINARY VARIABLES
538	!
539
540	if (IX .lt. im) stop
541	!
542	IPRT = 0
543	IF(IPRT.EQ.1) THEN
544	!!! LATD = 0
545	LOND = 0
546	ELSE
547	!!! LATD = 0
548	LOND = 0
549	ENDIF
550	!
551	gravi = 1.0 / grav
552	DT = 2. * DELTIM
553	RDT = 1. / DT
554	KMPBL = KM / 2
555	!
556	do k=1,km
557	do i=1,im
558	zi(i,k) = phii(i,k) * gravi
559	zl(i,k) = phil(i,k) * gravi
560	enddo
561	enddo
562	!
563	do k=1,kmpbl
564	do i=1,im
565	theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
566	enddo
567	enddo
568	!
569	DO K = 1,KM-1
570	DO I=1,IM
571	RDZT(I,K) = GOR * PRSI(I,K+1) / (PRSL(I,K) - PRSL(I,K+1))
572	ENDDO
573	ENDDO
574	!
575	DO I = 1,IM
576	DUSFC(I) = 0.
577	DVSFC(I) = 0.
578	DTSFC(I) = 0.
579	DQSFC(I) = 0.
580	HGAMT(I) = 0.
581	HGAMQ(I) = 0.
582	WSCALE(I) = 0.
583	KPBL(I) = 1
584	HPBL(I) = ZI(I,2)
585	PBLFLG(I) = .TRUE.
586	SFCFLG(I) = .TRUE.
587	IF(RBSOIL(I).GT.0.0) SFCFLG(I) = .FALSE.
588	ENDDO
589	!!
590	DO I=1,IM
591	RDZT1 = GOR * prSL(i,1) / DEL(i,1)
592	! BET1 = DTRDZT1SPD1(I)/T1(I,1)
593	BETA(I) = DT*RDZT1/T1(I,1)
594	! BETAW(I) = BET1*CD(I)
595	! BETAT(I) = BET1*CH(I)
596	! BETAQ(I) = DPHI(I)*BETAT(I)
597	ENDDO
598	!
599	DO I=1,IM
600	! ZL1(i) = 0.-(T1(I,1)+TSEA(I))/2.LOG(PRSL(I,1)/PRSI(I,1))ROG
601	! USTAR(I) = SQRT(CD(I)SPD1(I)*2)
602	USTAR(I) = SQRT(STRESS(I))
603	ENDDO
604	!
605	DO I=1,IM
606	THESV(I) = TSEA(I)(1.+FVMAX(QSS(I),QMIN))
607	THE1(I) = THETA(I,1)
608	THE1V(I) = THE1(I)(1.+FVMAX(Q1(I,1,1),QMIN))
609	THERMAL(I) = THE1V(I)
610	! DTHE1 = (THE1(I)-TSEA(I))
611	! DQ1 = (MAX(Q1(I,1,1),QMIN) - MAX(QSS(I),QMIN))
612	! HEAT(I) = -CH(I)SPD1(I)DTHE1
613	! EVAP(I) = -CH(I)SPD1(I)DQ1
614	ENDDO
615	!
616	!
617	! COMPUTE THE FIRST GUESS OF PBL HEIGHT
618	!
619	DO I=1,IM
620	STABLE(I) = .FALSE.
621	! ZL(i,1) = ZL1(i)
622	RBUP(I) = RBSOIL(I)
623	ENDDO
624	DO K = 2, KMPBL
625	DO I = 1, IM
626	IF(.NOT.STABLE(I)) THEN
627	RBDN(I) = RBUP(I)
628	! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
629	! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
630	THEKV(I) = THETA(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
631	SPDK2 = MAX(RCL(i)(U1(i,k)U1(i,k)+V1(i,k)*V1(i,k)),RONE)
632	RBUP(I) = (THEKV(I)-THE1V(I))(GZL(I,k)/THE1V(I))/SPDK2
633	KPBL(I) = K
634	STABLE(I) = RBUP(I).GT.RBCR
635	ENDIF
636	ENDDO
637	ENDDO
638	!
639	DO I = 1,IM
640	K = KPBL(I)
641	IF(RBDN(I).GE.RBCR) THEN
642	RBINT = 0.
643	ELSEIF(RBUP(I).LE.RBCR) THEN
644	RBINT = 1.
645	ELSE
646	RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I))
647	ENDIF
648	HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,K)-ZL(I,K-1))
649	IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
650	ENDDO
651	!!
652	DO I=1,IM
653	HOL = MAX(RBSOIL(I)FM(I)FM(I)/FH(I),RIMIN)
654	IF(SFCFLG(I)) THEN
655	HOL = MIN(HOL,-ZFMIN)
656	ELSE
657	HOL = MAX(HOL,ZFMIN)
658	ENDIF
659	!
660	! HOL = HOLHPBL(I)/ZL1(I)SFCFRAC
661	HOL = HOLHPBL(I)/ZL(I,1)SFCFRAC
662	IF(SFCFLG(I)) THEN
663	! PHIM = (1.-APHI16HOL)*(-1./4.)
664	! PHIH = (1.-APHI16HOL)*(-1./2.)
665	TEM = 1.0 / (1. - APHI16*HOL)
666	PHIH(I) = SQRT(TEM)
667	PHIM(I) = SQRT(PHIH(I))
668	ELSE
669	PHIM(I) = (1.+APHI5*HOL)
670	PHIH(I) = PHIM(I)
671	ENDIF
672	WSCALE(I) = USTAR(I)/PHIM(I)
673	WSCALE(I) = MIN(WSCALE(I),USTAR(I)*APHI16)
674	WSCALE(I) = MAX(WSCALE(I),USTAR(I)/APHI5)
675	ENDDO
676	!
677	! COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
678	! UNDER UNSTABLE CONDITIONS
679	!
680	DO I = 1,IM
681	SFLUX = HEAT(I) + EVAP(I)FVTHE1(I)
682	IF(SFCFLG(I).AND.SFLUX.GT.0.0) THEN
683	HGAMT(I) = MIN(CFAC*HEAT(I)/WSCALE(I),GAMCRT)
684	HGAMQ(I) = MIN(CFAC*EVAP(I)/WSCALE(I),GAMCRQ)
685	VPERT = HGAMT(I) + FVTHE1(I)HGAMQ(I)
686	VPERT = MIN(VPERT,GAMCRT)
687	THERMAL(I) = THERMAL(I) + MAX(VPERT,RZERO)
688	HGAMT(I) = MAX(HGAMT(I),RZERO)
689	HGAMQ(I) = MAX(HGAMQ(I),RZERO)
690	ELSE
691	PBLFLG(I) = .FALSE.
692	ENDIF
693	ENDDO
694	!
695	DO I = 1,IM
696	IF(PBLFLG(I)) THEN
697	KPBL(I) = 1
698	HPBL(I) = ZI(I,2)
699	ENDIF
700	ENDDO
701	!
702	! ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
703	!
704	DO I = 1, IM
705	IF(PBLFLG(I)) THEN
706	STABLE(I) = .FALSE.
707	RBUP(I) = RBSOIL(I)
708	ENDIF
709	ENDDO
710	DO K = 2, KMPBL
711	DO I = 1, IM
712	IF(.NOT.STABLE(I).AND.PBLFLG(I)) THEN
713	RBDN(I) = RBUP(I)
714	! ZL(I,k) = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
715	! & LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
716	THEKV(I) = THETA(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
717	SPDK2 = MAX(RCL(i)(U1(i,k)U1(i,k)+V1(i,k)*V1(i,k)),RONE)
718	RBUP(I) = (THEKV(I)-THERMAL(I))(GZL(I,k)/THE1V(I))/SPDK2
719	KPBL(I) = K
720	STABLE(I) = RBUP(I).GT.RBCR
721	ENDIF
722	ENDDO
723	ENDDO
724	!
725	DO I = 1,IM
726	IF(PBLFLG(I)) THEN
727	K = KPBL(I)
728	IF(RBDN(I).GE.RBCR) THEN
729	RBINT = 0.
730	ELSEIF(RBUP(I).LE.RBCR) THEN
731	RBINT = 1.
732	ELSE
733	RBINT = (RBCR-RBDN(I))/(RBUP(I)-RBDN(I))
734	ENDIF
735	HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,k)-ZL(I,K-1))
736	IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
737	IF(KPBL(I).LE.1) PBLFLG(I) = .FALSE.
738	ENDIF
739	ENDDO
740	!!
741	!
742	! COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
743	!
744	DO K = 1, KMPBL
745	DO I=1,IM
746	IF(KPBL(I).GT.K) THEN
747	PRINV = 1.0 / (PHIH(I)/PHIM(I)+CFACVK.1)
748	PRINV = MIN(PRINV,PRMAX)
749	PRINV = MAX(PRINV,PRMIN)
750	! ZFAC = MAX((1.-(ZI(I,K+1)-ZL1(I))/ &
751	! & (HPBL(I)-ZL1(I))), ZFMIN)
752	ZFAC = MAX((1.-(ZI(I,K+1)-ZL(I,1))/ &
753	& (HPBL(I)-ZL(I,1))), ZFMIN)
754	DKU(i,k) = XKZO + WSCALE(I)VKZI(I,K+1) &
755	& * ZFAC**PFAC
756	DKT(i,k) = DKU(i,k)*PRINV
757	DKO(i,k) = (DKU(i,k)-XKZO)*PRINV
758	DKU(i,k) = MIN(DKU(i,k),DKMAX)
759	DKU(i,k) = MAX(DKU(i,k),DKMIN)
760	DKT(i,k) = MIN(DKT(i,k),DKMAX)
761	DKT(i,k) = MAX(DKT(i,k),DKMIN)
762	DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
763	ENDIF
764	ENDDO
765	ENDDO
766	!
767	! COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
768	!
769	DO K = 1, KM-1
770	DO I=1,IM
771	IF(K.GE.KPBL(I)) THEN
772	! TI = 0.5*(T1(i,k)+T1(i,K+1))
773	TI = 2.0 / (T1(i,k)+T1(i,K+1))
774	! RDZ = RDZT(I,K)/TI
775	RDZ = RDZT(I,K) * TI
776	! RDZ = RDZT(I,K)
777	DW2 = RCL(i)((U1(i,k)-U1(i,K+1))*2 &
778	& + (V1(i,k)-V1(i,K+1))**2)
779	SHR2 = MAX(DW2,DW2MIN)RDZ*2
780	TVD = T1(i,k)(1.+FVMAX(Q1(i,k,1),QMIN))
781	TVU = T1(i,K+1)(1.+FVMAX(Q1(i,K+1,1),QMIN))
782	! BVF2 = G(GOCP+RDZ(TVU-TVD))/TI
783	BVF2 = G(GOCP+RDZ(TVU-TVD)) * TI
784	RI = MAX(BVF2/SHR2,RIMIN)
785	ZK = VK*ZI(I,K+1)
786	! RL2 = (ZKRLAM/(RLAM+ZK))*2
787	! DK = RL2*SQRT(SHR2)
788	RL2 = ZK*RLAM/(RLAM+ZK)
789	DK = RL2RL2SQRT(SHR2)
790	IF(RI.LT.0.) THEN ! UNSTABLE REGIME
791	SRI = SQRT(-RI)
792	DKU(i,k) = XKZO + DK(1+8.(-RI)/(1+1.746*SRI))
793	! DKT(i,k) = XKZO + DK(1+8.(-RI)/(1+1.286*SRI))
794	tem = DK(1+8.(-RI)/(1+1.286*SRI))
795	DKT(i,k) = XKZO + tem
796	DKO(i,k) = tem
797	ELSE ! STABLE REGIME
798	! DKT(i,k) = XKZO + DK/(1+5.RI)*2
799	tem = DK/(1+5.RI)*2
800	DKT(i,k) = XKZO + tem
801	DKO(i,k) = tem
802	PRNUM = 1.0 + 2.1*RI
803	PRNUM = MIN(PRNUM,PRMAX)
804	DKU(i,k) = (DKT(i,k)-XKZO)*PRNUM + XKZO
805	ENDIF
806	!
807	DKU(i,k) = MIN(DKU(i,k),DKMAX)
808	DKU(i,k) = MAX(DKU(i,k),DKMIN)
809	DKT(i,k) = MIN(DKT(i,k),DKMAX)
810	DKT(i,k) = MAX(DKT(i,k),DKMIN)
811	DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
812	!
813	!!! IF(I.EQ.LOND.AND.LAT.EQ.LATD) THEN
814	!!! PRNUM = DKU(k)/DKT(k)
815	!!! WRITE(IUN,610) K,PRNUM,DKT(k),DKU(k),RL2,RI,
816	!!! 1 BVF2,SHR2
817	!!! ENDIF
818	!
819	ENDIF
820	ENDDO
821	ENDDO
822	!
823	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
824	!
825	DO I=1,IM
826	AD(I,1) = 1.
827	A1(I,1) = T1(i,1) + BETA(i) * HEAT(I)
828	A2(I,1) = Q1(i,1,1) + BETA(i) * EVAP(I)
829	! A1(I,1) = T1(i,1)-BETAT(I)*(THETA(i,1)-TSEA(I))
830	! A2(I,1) = Q1(i,1,1)-BETAQ(I)*
831	! & (MAX(Q1(i,1,1),QMIN)-MAX(QSS(I),QMIN))
832	ENDDO
833	!
834	DO K = 1,KM-1
835	DO I = 1,IM
836	DTODSD = DT/DEL(I,K)
837	DTODSU = DT/DEL(I,K+1)
838	DSIG = PRSL(I,K)-PRSL(I,K+1)
839	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
840	! RDZ = RDZT(I,K)
841	tem1 = DSIG * DKT(i,k) * RDZ
842	IF(PBLFLG(I).AND.K.LT.KPBL(I)) THEN
843	! DSDZT = DSIGDKT(i,k)RDZ*(GOCP-HGAMT(I)/HPBL(I))
844	! DSDZQ = DSIGDKT(i,k)RDZ*(-HGAMQ(I)/HPBL(I))
845	tem = 1.0 / HPBL(I)
846	DSDZT = tem1 * (GOCP-HGAMT(I)*tem)
847	DSDZQ = tem1 * (-HGAMQ(I)*tem)
848	A2(I,k) = A2(I,k)+DTODSD*DSDZQ
849	A2(I,k+1) = Q1(i,k+1,1)-DTODSU*DSDZQ
850	ELSE
851	! DSDZT = DSIGDKT(i,k)RDZ*(GOCP)
852	DSDZT = tem1 * GOCP
853	A2(I,k+1) = Q1(i,k+1,1)
854	ENDIF
855	! DSDZ2 = DSIGDKT(i,k)RDZ*RDZ
856	DSDZ2 = tem1 * RDZ
857	AU(I,k) = -DTODSD*DSDZ2
858	AL(I,k) = -DTODSU*DSDZ2
859	AD(I,k) = AD(I,k)-AU(I,k)
860	AD(I,k+1) = 1.-AL(I,k)
861	A1(I,k) = A1(I,k)+DTODSD*DSDZT
862	A1(I,k+1) = T1(i,k+1)-DTODSU*DSDZT
863	ENDDO
864	ENDDO
865	!
866	! SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
867	!
868	CALL TRIDIN(IM,KM,1,AL,AD,AU,A1,A2,AU,A1,A2)
869	!
870	! RECOVER TENDENCIES OF HEAT AND MOISTURE
871	!
872	DO K = 1,KM
873	DO I = 1,IM
874	TTEND = (A1(I,k)-T1(i,k))*RDT
875	QTEND = (A2(I,k)-Q1(i,k,1))*RDT
876	TAU(i,k) = TAU(i,k)+TTEND
877	RTG(I,k,1) = RTG(i,k,1)+QTEND
878	DTSFC(I) = DTSFC(I)+CONTDEL(I,K)TTEND
879	DQSFC(I) = DQSFC(I)+CONQDEL(I,K)QTEND
880	ENDDO
881	ENDDO
882	!
883	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
884	!
885	DO I=1,IM
886	! AD(I,1) = 1.+BETAW(I)
887	AD(I,1) = 1.0 + BETA(i) * STRESS(I) / SPD1(I)
888	A1(I,1) = U1(i,1)
889	A2(I,1) = V1(i,1)
890	! AD(I,1) = 1.0
891	! tem = 1.0 + BETA(I) * STRESS(I) / SPD1(I)
892	! A1(I,1) = U1(i,1) * tem
893	! A2(I,1) = V1(i,1) * tem
894	ENDDO
895	!
896	DO K = 1,KM-1
897	DO I=1,IM
898	DTODSD = DT/DEL(I,K)
899	DTODSU = DT/DEL(I,K+1)
900	DSIG = PRSL(I,K)-PRSL(I,K+1)
901	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,k+1))
902	! RDZ = RDZT(I,K)
903	DSDZ2 = DSIGDKU(i,k)RDZ*RDZ
904	AU(I,k) = -DTODSD*DSDZ2
905	AL(I,k) = -DTODSU*DSDZ2
906	AD(I,k) = AD(I,k)-AU(I,k)
907	AD(I,k+1) = 1.-AL(I,k)
908	A1(I,k+1) = U1(i,k+1)
909	A2(I,k+1) = V1(i,k+1)
910	ENDDO
911	ENDDO
912	!
913	! SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
914	!
915	CALL TRIDI2(IM,KM,AL,AD,AU,A1,A2,AU,A1,A2)
916	!
917	! RECOVER TENDENCIES OF MOMENTUM
918	!
919	DO K = 1,KM
920	DO I = 1,IM
921	CONWRC = CONW*SQRT(RCL(i))
922	UTEND = (A1(I,k)-U1(i,k))*RDT
923	VTEND = (A2(I,k)-V1(i,k))*RDT
924	DU(i,k) = DU(i,k)+UTEND
925	DV(i,k) = DV(i,k)+VTEND
926	DUSFC(I) = DUSFC(I)+CONWRCDEL(I,K)UTEND
927	DVSFC(I) = DVSFC(I)+CONWRCDEL(I,K)VTEND
928	ENDDO
929	ENDDO
930	!!
931	!
932	! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR TRACERS
933	!
934	if (ntrac .ge. 2) then
935	DO I=1,IM
936	AD(I,1) = 1.
937	ENDDO
938	do k = 2, ntrac
939	is = (k-2) * km
940	do i = 1, im
941	AT(I,1+is) = Q1(i,1,k)
942	enddo
943	enddo
944	!
945	DO K = 1,KM-1
946	DO I = 1,IM
947	DTODSD = DT/DEL(I,K)
948	DTODSU = DT/DEL(I,K+1)
949	DSIG = PRSL(I,K)-PRSL(I,K+1)
950	RDZ = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
951	tem1 = DSIG * DKT(i,k) * RDZ
952	DSDZ2 = tem1 * RDZ
953	AU(I,k) = -DTODSD*DSDZ2
954	AL(I,k) = -DTODSU*DSDZ2
955	AD(I,k) = AD(I,k)-AU(I,k)
956	AD(I,k+1) = 1.-AL(I,k)
957	ENDDO
958	ENDDO
959	do kk = 2, ntrac
960	is = (kk-2) * km
961	do k = 1, km - 1
962	do i = 1, im
963	AT(I,k+1+is) = Q1(i,k+1,kk)
964	enddo
965	enddo
966	enddo
967	!
968	! SOLVE TRIDIAGONAL PROBLEM FOR TRACERS
969	!
970	CALL TRIDIT(IM,KM,ntrac-1,AL,AD,AU,AT,AU,AT)
971	!
972	! RECOVER TENDENCIES OF TRACERS
973	!
974	do kk = 2, ntrac
975	is = (kk-2) * km
976	do k = 1, km
977	do i = 1, im
978	QTEND = (AT(I,K+is)-Q1(i,K,kk))*RDT
979	RTG(i,K,kk) = RTG(i,K,kk) + QTEND
980	enddo
981	enddo
982	enddo
983	endif
984	!!
985	RETURN
986	END SUBROUTINE MONINP
987	!FPP$ NOCONCUR R
988	!-----------------------------------------------------------------------
989	SUBROUTINE TRIDI2(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
990	!sela %INCLUDE DBTRIDI2;
991	!
992	USE MODULE_GFS_MACHINE, ONLY : kind_phys
993	implicit none
994	integer k,n,l,i
995	real(kind=kind_phys) fk
996	!
997	real(kind=kind_phys) CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
998	& AU(L,N-1),A1(L,N),A2(L,N)
999	!-----------------------------------------------------------------------
1000	DO I=1,L
1001	FK = 1./CM(I,1)
1002	AU(I,1) = FK*CU(I,1)
1003	A1(I,1) = FK*R1(I,1)
1004	A2(I,1) = FK*R2(I,1)
1005	ENDDO
1006	DO K=2,N-1
1007	DO I=1,L
1008	FK = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1009	AU(I,K) = FK*CU(I,K)
1010	A1(I,K) = FK(R1(I,K)-CL(I,K)A1(I,K-1))
1011	A2(I,K) = FK(R2(I,K)-CL(I,K)A2(I,K-1))
1012	ENDDO
1013	ENDDO
1014	DO I=1,L
1015	FK = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1016	A1(I,N) = FK(R1(I,N)-CL(I,N)A1(I,N-1))
1017	A2(I,N) = FK(R2(I,N)-CL(I,N)A2(I,N-1))
1018	ENDDO
1019	DO K=N-1,1,-1
1020	DO I=1,L
1021	A1(I,K) = A1(I,K)-AU(I,K)*A1(I,K+1)
1022	A2(I,K) = A2(I,K)-AU(I,K)*A2(I,K+1)
1023	ENDDO
1024	ENDDO
1025	!-----------------------------------------------------------------------
1026	RETURN
1027	END SUBROUTINE TRIDI2
1028	!FPP$ NOCONCUR R
1029	!-----------------------------------------------------------------------
1030	SUBROUTINE TRIDIN(L,N,nt,CL,CM,CU,R1,R2,AU,A1,A2)
1031	!sela %INCLUDE DBTRIDI2;
1032	!
1033	USE MODULE_GFS_MACHINE, ONLY : kind_phys
1034	implicit none
1035	integer is,k,kk,n,nt,l,i
1036	real(kind=kind_phys) fk(L)
1037	!
1038	real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1), &
1039	& R1(L,N), R2(L,N*nt), &
1040	& AU(L,N-1), A1(L,N), A2(L,N*nt), &
1041	& FKK(L,2:N-1)
1042	!-----------------------------------------------------------------------
1043	DO I=1,L
1044	FK(I) = 1./CM(I,1)
1045	AU(I,1) = FK(I)*CU(I,1)
1046	A1(I,1) = FK(I)*R1(I,1)
1047	ENDDO
1048	do k = 1, nt
1049	is = (k-1) * n
1050	do i = 1, l
1051	a2(i,1+is) = fk(I) * r2(i,1+is)
1052	enddo
1053	enddo
1054	DO K=2,N-1
1055	DO I=1,L
1056	FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1057	AU(I,K) = FKK(I,K)*CU(I,K)
1058	A1(I,K) = FKK(I,K)(R1(I,K)-CL(I,K)A1(I,K-1))
1059	ENDDO
1060	ENDDO
1061	do kk = 1, nt
1062	is = (kk-1) * n
1063	DO K=2,N-1
1064	DO I=1,L
1065	A2(I,K+is) = FKK(I,K)(R2(I,K+is)-CL(I,K)A2(I,K+is-1))
1066	ENDDO
1067	ENDDO
1068	ENDDO
1069	DO I=1,L
1070	FK(I) = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1071	A1(I,N) = FK(I)(R1(I,N)-CL(I,N)A1(I,N-1))
1072	ENDDO
1073	do k = 1, nt
1074	is = (k-1) * n
1075	do i = 1, l
1076	A2(I,N+is) = FK(I)(R2(I,N+is)-CL(I,N)A2(I,N+is-1))
1077	enddo
1078	enddo
1079	DO K=N-1,1,-1
1080	DO I=1,L
1081	A1(I,K) = A1(I,K) - AU(I,K)*A1(I,K+1)
1082	ENDDO
1083	ENDDO
1084	do kk = 1, nt
1085	is = (kk-1) * n
1086	DO K=n-1,1,-1
1087	DO I=1,L
1088	A2(I,K+is) = A2(I,K+is) - AU(I,K)*A2(I,K+is+1)
1089	ENDDO
1090	ENDDO
1091	ENDDO
1092	!-----------------------------------------------------------------------
1093	RETURN
1094	END SUBROUTINE TRIDIN
1095	SUBROUTINE TRIDIT(L,N,nt,CL,CM,CU,RT,AU,AT)
1096	!sela %INCLUDE DBTRIDI2;
1097	!
1098	USE MODULE_GFS_MACHINE, ONLY : kind_phys
1099	implicit none
1100	integer is,k,kk,n,nt,l,i
1101	real(kind=kind_phys) fk(L)
1102	!
1103	real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1), &
1104	& RT(L,N*nt), &
1105	& AU(L,N-1), AT(L,N*nt), &
1106	& FKK(L,2:N-1)
1107	!-----------------------------------------------------------------------
1108	DO I=1,L
1109	FK(I) = 1./CM(I,1)
1110	AU(I,1) = FK(I)*CU(I,1)
1111	ENDDO
1112	do k = 1, nt
1113	is = (k-1) * n
1114	do i = 1, l
1115	at(i,1+is) = fk(I) * rt(i,1+is)
1116	enddo
1117	enddo
1118	DO K=2,N-1
1119	DO I=1,L
1120	FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1121	AU(I,K) = FKK(I,K)*CU(I,K)
1122	ENDDO
1123	ENDDO
1124	do kk = 1, nt
1125	is = (kk-1) * n
1126	DO K=2,N-1
1127	DO I=1,L
1128	AT(I,K+is) = FKK(I,K)(RT(I,K+is)-CL(I,K)AT(I,K+is-1))
1129	ENDDO
1130	ENDDO
1131	ENDDO
1132	DO I=1,L
1133	FK(I) = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1134	ENDDO
1135	do k = 1, nt
1136	is = (k-1) * n
1137	do i = 1, l
1138	AT(I,N+is) = FK(I)(RT(I,N+is)-CL(I,N)AT(I,N+is-1))
1139	enddo
1140	enddo
1141	do kk = 1, nt
1142	is = (kk-1) * n
1143	DO K=n-1,1,-1
1144	DO I=1,L
1145	AT(I,K+is) = AT(I,K+is) - AU(I,K)*AT(I,K+is+1)
1146	ENDDO
1147	ENDDO
1148	ENDDO
1149	!-----------------------------------------------------------------------
1150	RETURN
1151	END SUBROUTINE TRIDIT
1152
1153	!-----------------------------------------------------------------------
1154
1155	END MODULE module_bl_gfs

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

Download in other formats:

Original Format