PageRenderTime 65ms CodeModel.GetById 28ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_bl_gfs.F

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