Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_gfs.F @ 3567

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

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

Download in other formats:

Original Format