| blob | f747205c0af011c3108df430e38469171cb19afd |
1 !LWRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_bl_gfs
5 CONTAINS
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 !-------------------------------------------------------------------
91 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
92 ims,ime, jms,jme, kms,kme, &
93 its,ite, jts,jte, kts,kte
95 LOGICAL, INTENT(IN) :: flag_qi
97 REAL, INTENT(IN) :: &
98 CP, &
99 DT, &
100 EP1, &
101 G, &
102 KARMAN, &
103 R, &
104 ROVCP, &
105 ROVG
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
120 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
121 RTHBLTEN, &
122 RQCBLTEN, &
123 RQVBLTEN, &
124 RUBLTEN, &
125 RVBLTEN
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
137 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
138 PBL, &
139 UST, &
140 WSPD
142 INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
143 KPBL2D
145 !--------------------------- OPTIONAL VARS ------------------------------
146 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
147 OPTIONAL, INTENT(IN) :: &
148 QI3D
150 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
151 OPTIONAL, INTENT(INOUT) :: &
152 RQIBLTEN
154 !--------------------------- LOCAL VARS ------------------------------
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
169 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
170 PHII, &
171 PRSI
173 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 3) :: &
174 Q1, &
175 RTG
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
197 REAL (kind=kind_phys) :: &
198 CPM, &
199 DELTIM, &
200 FMTMP, &
201 RRHOX
203 INTEGER, DIMENSION( its:ite ) :: &
204 KPBL
206 INTEGER :: &
207 I, &
208 IM, &
209 J, &
210 K, &
211 KM, &
212 KTEM, &
213 KTEP, &
214 KX, &
215 L, &
216 NTRAC
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
227 DO J=jts,jte
229 DO i=its,ite
230 RRHOX=(R*T3D(I,KTS,J)*(1.+EP1*QV3D(I,KTS,J)))/PSFC(I,J)
231 CPM=CP*(1.+0.8*QV3D(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)=KARMAN*KARMAN*WSPD(i,j)*WSPD(i,j)/(FMTMP*FMTMP)
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
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
262 DO k=kts,kte
263 DO i=its,ite
264 PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
265 ENDDO
266 ENDDO
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
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
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
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
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)
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
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
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
333 ENDDO
336 END SUBROUTINE BL_GFS
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
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
364 jtf=min0(jte,jde-1)
365 ktf=min0(kte,kde-1)
366 itf=min0(ite,ide-1)
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
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
392 END SUBROUTINE gfsinit
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
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 !
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 = DT*RDZT1*SPD1(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.+FV*MAX(QSS(I),QMIN))
574 THE1(I) = THETA(I,1)
575 THE1V(I) = THE1(I)*(1.+FV*MAX(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.+FV*MAX(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))*(G*ZL(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 = HOL*HPBL(I)/ZL1(I)*SFCFRAC
628 HOL = HOL*HPBL(I)/ZL(I,1)*SFCFRAC
629 IF(SFCFLG(I)) THEN
630 ! PHIM = (1.-APHI16*HOL)**(-1./4.)
631 ! PHIH = (1.-APHI16*HOL)**(-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)*FV*THE1(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) + FV*THE1(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.+FV*MAX(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))*(G*ZL(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)+CFAC*VK*.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)*VK*ZI(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.+FV*MAX(Q1(i,k,1),QMIN))
748 TVU = T1(i,K+1)*(1.+FV*MAX(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 = (ZK*RLAM/(RLAM+ZK))**2
754 ! DK = RL2*SQRT(SHR2)
755 RL2 = ZK*RLAM/(RLAM+ZK)
756 DK = RL2*RL2*SQRT(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 = DSIG*DKT(i,k)*RDZ*(GOCP-HGAMT(I)/HPBL(I))
811 ! DSDZQ = DSIG*DKT(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 = DSIG*DKT(i,k)*RDZ*(GOCP)
819 DSDZT = tem1 * GOCP
820 A2(I,k+1) = Q1(i,k+1,1)
821 ENDIF
822 ! DSDZ2 = DSIG*DKT(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)+CONT*DEL(I,K)*TTEND
846 DQSFC(I) = DQSFC(I)+CONQ*DEL(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 = DSIG*DKU(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)+CONWRC*DEL(I,K)*UTEND
894 DVSFC(I) = DVSFC(I)+CONWRC*DEL(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
1120 !-----------------------------------------------------------------------
1122 END MODULE module_bl_gfs