CMake netCDF Compatibility with WPS (#2121)
[WRF.git] / phys / module_bl_eepsilon.F
blob58b5f7b7393dad6ea7cc05d1d89456562c2040a2
1 !WRF:model_layer:physics\r
2 !\r
3 !####################e-epsilon scheme##################################\r
4 !   taken from the IPRC IRAM - Yuqing Wang, university of hawaii\r
5 !   added and improved by Chunxi Zhang and Yuqing Wang to wrf4.2, 2020\r
6 !   refenrences: Langland and Liou (1996, mwr, 124, 905-918)\r
7 !                Zhang et al.      (2020, mwr, 148, 1121-1145)\r
8 !######################################################################\r
9 \r
10 !  Updates since the original release\r
11 !  03/2023\r
12 !  Added the cloud top cooling effect for stratocumulus, which is similar\r
13 !      to the one used in the MYNN PBL scheme. The critial step is to\r
14 !      idendity stratocumulus. We used the method by Marquet and Bechtold (2020)\r
15 !  Added addtinal source term for stable boundary layer in the epsilon \r
16 !      equation to represent the energy drain on the eddy scale\r
17 !      (Zeng et al. (2020))\r
18 !  Fixed a bug for the TKE global varaible. A new variable pek_pbl is defined \r
19 !  instead of  using tke_pbl in Register.EM_COMMON. Now pek_pbl,pep_pbl,pek_adv\r
20 !  and pep_adv are all on the half levels. \r
21 !  Addded a few tweaks\r
22 !\r
23 !######################################################################\r
24 module module_bl_eeps\r
26    use module_model_constants, only:rgas=>r_d,rv=>r_v, &\r
27    &       t13=>Prandtl,ep1=>EP_1,ep2=>EP_2,grv=>g, &\r
28    &       akv=>KARMAN,cp,rcp,xlv,xls,xlf,p1000mb\r
29      \r
30    real, private, parameter :: cs2=17.67,cs3=273.15,cs4=29.65\r
31    real, private, parameter :: ci2=22.514,ci3=273.16,ci4=0.0\r
32    real, private, parameter :: hgfr=233.15,t0=273.15\r
33    real, private, parameter :: epsl=1.e-20\r
34    real, private, parameter :: us2min=0.1\r
35    real, private, parameter :: akvmx = 1600.\r
36    real, private, parameter :: ricr=0.25,bt=5.0\r
37    real, private, parameter :: minpek = 1.0e-4\r
38    real, private, parameter :: minpep = 1.0e-7\r
39    real, private, parameter :: maxpek = 160.\r
40    real, private, parameter :: maxpep = 25.0\r
41    real, private, parameter :: zfmin = 0.01\r
42    integer,private,parameter:: bl_eeps_tkeadv = 1\r
43 !Adding top-down diffusion driven by cloud-top radiative cooling\r
44    integer,private,parameter:: bl_eeps_topdown = 1\r
45 !---------------------------------------------------------------\r
46 ! ci--- five constants used in TKE tubrbulence closure scheme\r
47 !---------------------------------------------------------------\r
48 ! Detering and Etling (1986); Langland and Liou (1996)\r
49 !   real,parameter :: c1=1.35,c2=0.026,c3=1.13,c4=1.90,c5=0.77\r
50 ! Duynkerke and Driedonks (1987); Stubley and Rooney (1986)\r
51    real,parameter :: c1=1.35,c2=0.09,c3=1.44,c4=1.92,c5=0.77\r
52 ! Beljaars et al. (1987)\r
53 !   real,parameter :: c1=1.35,c2=0.032,c3=1.44,c4=1.92,c5=0.54\r
54 !\r
55 contains\r
56 !-------------------------------------------------------------------------------\r
57 !\r
58    subroutine eeps(u3d,v3d,t3d,qv3d,qc3d,qi3d,qr3d,qs3d,qg3d,p3d,pi3d,    &\r
59                   rho3d,rthraten,rublten,rvblten,rthblten,                &\r
60                   rqvblten,rqcblten,rqiblten,                             &\r
61                   pek,pep,                                                &\r
62                   dz8w,psfc,qsfc,tsk,ust,rmol,wspd,                       &\r
63                   xland,hfx,qfx,                                          &\r
64                   dt,dx,itimestep,exch_h,exch_m,pblh,kpbl,                &\r
65                   pek_adv,pep_adv,                                        &\r
66                   ids,ide, jds,jde, kds,kde,                              &\r
67                   ims,ime, jms,jme, kms,kme,                              &\r
68                   its,ite, jts,jte, kts,kte                               &\r
69                                                       )\r
70 !-------------------------------------------------------------------------------\r
71    implicit none\r
72 !-------------------------------------------------------------------------------\r
73 !-- u3d         3d u-velocity interpolated to theta points (m/s)\r
74 !-- v3d         3d v-velocity interpolated to theta points (m/s)\r
75 !-- t3d         temperature (k)\r
76 !-- qv3d        3d water vapor mixing ratio (kg/kg)\r
77 !-- qc3d        3d cloud mixing ratio (kg/kg)\r
78 !-- qi3d        3d ice mixing ratio (kg/kg)\r
79 !-- qr3d        3d rain mixing ratio (kg/kg)\r
80 !-- qs3d        3d snow mixing ratio (kg/kg)\r
81 !-- qg3d        3d graupel mixing ratio (kg/kg)\r
82 !               (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)\r
83 !-- p3d         3d pressure (pa)\r
84 !-- pi3d        3d exner function (dimensionless)\r
85 !-- rho3d       3d air density (kg/m^3)\r
86 !-- rthraten    theta tendency due to\r
87 !               radiation parameterization (K/s)\r
88 !-- rublten     u tendency due to\r
89 !               pbl parameterization (m/s/s)\r
90 !-- rvblten     v tendency due to\r
91 !               pbl parameterization (m/s/s)\r
92 !-- rthblten    theta tendency due to\r
93 !               pbl parameterization (K/s)\r
94 !-- rqvblten    qv tendency due to\r
95 !               pbl parameterization (kg/kg/s)\r
96 !-- rqcblten    qc tendency due to\r
97 !               pbl parameterization (kg/kg/s)\r
98 !-- rqiblten    qi tendency due to\r
99 !               pbl parameterization (kg/kg/s)\r
100 !-- dz8w        dz between full levels (m)\r
101 !-- pek         3d turbulent kinetic energy (TKE, m2/s2)\r
102 !-- pep         3d dissipation rate of TKE \r
103 !-- psfc        pressure at the surface (pa)\r
104 !-- qsfc        ground saturated mixing ratio\r
105 !-- ust         u* in similarity theory (m/s)\r
106 !-- rmol        inverse (1/L) of Monin-Obukhov length (m)\r
107 !-- xland       land mask (1 for land, 2 for water)\r
108 !-- tsk         surface skin temperature\r
109 !-- hfx         upward heat flux at the surface (w/m^2)\r
110 !-- qfx         upward moisture flux at the surface (kg/m^2/s)\r
111 !-- dt          time step (s)\r
112 !-- ids         start index for i in domain\r
113 !-- ide         end index for i in domain\r
114 !-- jds         start index for j in domain\r
115 !-- jde         end index for j in domain\r
116 !-- kds         start index for k in domain\r
117 !-- kde         end index for k in domain\r
118 !-- ims         start index for i in memory\r
119 !-- ime         end index for i in memory\r
120 !-- jms         start index for j in memory\r
121 !-- jme         end index for j in memory\r
122 !-- kms         start index for k in memory\r
123 !-- kme         end index for k in memory\r
124 !-- its         start index for i in tile\r
125 !-- ite         end index for i in tile\r
126 !-- jts         start index for j in tile\r
127 !-- jte         end index for j in tile\r
128 !-- kts         start index for k in tile\r
129 !-- kte         end index for k in tile\r
130 !-------------------------------------------------------------------------------\r
132    integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde,                &\r
133                                      ims,ime, jms,jme, kms,kme,                &\r
134                                      its,ite, jts,jte, kts,kte\r
136    integer,  intent(in   )   ::      itimestep\r
138    real,     intent(in   )   ::      dt, dx\r
141    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
142              intent(in   )   ::                                          qv3d, &\r
143                                                                          qc3d, &\r
144                                                                          qi3d, &\r
145                                                                           p3d, &\r
146                                                                          pi3d, &\r
147                                                                           t3d, &\r
148                                                                          dz8w, &\r
149                                                                           u3d, &\r
150                                                                           v3d, &\r
151                                                                         rho3d\r
152    \r
153    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
154              optional                                                        , &\r
155              intent(in   )   ::                                          qr3d, &\r
156                                                                          qs3d, &\r
157                                                                          qg3d\r
159    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
160              intent(in)   ::                                         rthraten\r
162    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
163              intent(inout)   ::                                       rublten, &\r
164                                                                       rvblten, &\r
165                                                                      rthblten, &\r
166                                                                      rqvblten, &\r
167                                                                      rqcblten, &\r
168                                                                      rqiblten\r
170    real,     dimension( ims:ime, jms:jme )                                   , &\r
171              intent(in   )   ::                                         xland, &\r
172                                                                           hfx, &\r
173                                                                           qfx, &\r
174                                                                          psfc, &\r
175                                                                          qsfc, &\r
176                                                                           tsk, &\r
177                                                                           ust, &\r
178                                                                           rmol,&\r
179                                                                           wspd\r
181    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
182              intent(inout)   ::                                       pek, pep\r
184    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
185              intent(inout)   ::                                 pek_adv,pep_adv\r
187    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &\r
188              intent(out)   ::                                    exch_h,exch_m\r
190    real,     dimension( ims:ime, jms:jme )                                   , &\r
191              intent(out   )   ::                                          pblh\r
193    integer,  dimension( ims:ime, jms:jme )                                   , &\r
194              intent(out   )   ::                                          kpbl\r
196 !local\r
197       real,     dimension( its:ite, kts:kte+1)::              zi,ghti  \r
198       real,     dimension( its:ite, kts:kte ) ::              zl                                                   \r
199       \r
200       real,     dimension( its:ite, kts:kte ) ::              u1, &\r
201                                                               v1, &\r
202                                                               q1, &\r
203                                                               q2, &\r
204                                                               q3, &\r
205                                                               q4, &\r
206                                                               q5, &\r
207                                                               q6, &\r
208                                                               t1, &\r
209                                                               ghtl, &\r
210                                                               prsl, &\r
211                                                               rho1, &\r
212                                                          rthraten1\r
213       real,     dimension( its:ite, kts:kte ) ::          pek2d,pep2d\r
214       real,     dimension( its:ite, kts:kte+1 ) ::        exh2d,exm2d\r
215       real,     dimension( its:ite) ::               psfc1,ust1,rmol1,  &\r
216                                                      xland1,qsfc1,hfx1, &\r
217                                                      qfx1,tsk1,pblh1,wspd1\r
218       integer,  dimension( its:ite) ::               kpbl1\r
220       real    :: rdelt\r
221       integer ::  i,j,k,zz,im,kx,pp\r
223       im=ite-its+1\r
224       kx=kte-kts+1\r
225       rdelt = 1./dt\r
227       if(itimestep .eq. 1) then\r
228          pek_adv = pek\r
229          pep_adv = pep\r
230       end if\r
232       do j=jts,jte\r
233 ! --------------- compute zi and zl -----------------------------------------\r
234       do i=its,ite\r
235         zi(i,kts)=0.0\r
236       enddo\r
238       do k=kts,kte\r
239         do i=its,ite\r
240           zi(i,k+1)=zi(i,k)+dz8w(i,k,j)\r
241         enddo\r
242       enddo\r
244       do k=kts,kte\r
245         do i=its,ite\r
246           zl(i,k)=0.5*(zi(i,k)+zi(i,k+1))\r
247         enddo\r
248       enddo\r
250 ! reverse the vertical levels\r
251       pp = 0\r
252       do k=kts,kte\r
253         zz = kte-pp\r
254         do i=its,ite\r
255           u1(i,zz)=u3d(i,k,j)\r
256           v1(i,zz)=v3d(i,k,j)\r
257           t1(i,zz)=t3d(i,k,j)\r
258           q1(i,zz)=qv3d(i,k,j)/(1.+qv3d(i,k,j))\r
259           q2(i,zz)=qc3d(i,k,j)/(1.+qc3d(i,k,j))\r
260           q3(i,zz)=qi3d(i,k,j)/(1.+qi3d(i,k,j))\r
261           if(present(qr3d)) then\r
262             q4(i,zz)=qr3d(i,k,j)/(1.+qr3d(i,k,j))\r
263           else\r
264             q4(i,zz)=0.\r
265           end if\r
266           if(present(qs3d)) then\r
267             q5(i,zz)=qs3d(i,k,j)/(1.+qs3d(i,k,j))\r
268           else\r
269             q5(i,zz)=0.\r
270           end if\r
271           if(present(qg3d)) then\r
272             q6(i,zz)=qg3d(i,k,j)/(1.+qg3d(i,k,j))\r
273           else\r
274             q6(i,zz)=0.\r
275           end if\r
276           ghtl(i,zz)=zl(i,k)\r
277           prsl(i,zz) = p3d(i,k,j)\r
278           rho1(i,zz)=rho3d(i,k,j)\r
279           rthraten1(i,zz)=rthraten(i,k,j)\r
280           if(bl_eeps_tkeadv==1)then\r
281             pek2d(i,zz)= min(maxpek,max(minpek,pek_adv(i,k,j)))\r
282             pep2d(i,zz)= min(maxpep,max(minpep,pep_adv(i,k,j)))\r
283           else\r
284             pek2d(i,zz)= min(maxpek,max(minpek,pek(i,k,j)))\r
285             pep2d(i,zz)= min(maxpep,max(minpep,pep(i,k,j)))\r
286           end if\r
287         end do\r
288         pp = pp + 1\r
289       end do\r
291        pp = 0\r
292        do k=kts,kte+1\r
293         zz = kte+1-pp\r
294         do i=its,ite\r
295           ghti(i,zz) = zi(i,k)\r
296         enddo\r
297         pp = pp + 1\r
298       enddo\r
300       do i=its,ite\r
301         psfc1(i) = psfc(i,j)\r
302         ust1(i)  = ust(i,j)\r
303         rmol1(i) = rmol(i,j)\r
304         wspd1(i) = wspd(i,j)\r
305         xland1(i)= xland(i,j)\r
306         qsfc1(i) = qsfc(i,j)\r
307         hfx1(i)  = hfx(i,j)\r
308         qfx1(i)  = qfx(i,j)\r
309         tsk1(i)  = tsk(i,j)\r
310       end do\r
311 \r
312       call eeps2d(u1,v1,t1,q1,q2,q3,q4,q5,q6,prsl,ghtl,ghti,rho1,rthraten1,pek2d,pep2d &\r
313               ,psfc1,ust1,rmol1,wspd1,xland1,qsfc1,hfx1,qfx1,tsk1       &\r
314               ,dt,dx,itimestep,exh2d,exm2d,pblh1,kpbl1,im,kx)\r
316       pp = 0\r
317       do k=kts,kte\r
318         zz = kte-pp\r
319         do i=its,ite\r
320           rthblten(i,k,j)=(t1(i,zz)-t3d(i,k,j))/pi3d(i,k,j)*rdelt\r
321           rqvblten(i,k,j)=(q1(i,zz)/(1.-q1(i,zz))-qv3d(i,k,j))*rdelt\r
322           rqcblten(i,k,j)=(q2(i,zz)/(1.-q2(i,zz))-qc3d(i,k,j))*rdelt\r
323           rqiblten(i,k,j)=(q3(i,zz)/(1.-q3(i,zz))-qi3d(i,k,j))*rdelt\r
324           rublten(i,k,j) =(u1(i,zz)-u3d(i,k,j))*rdelt\r
325           rvblten(i,k,j) =(v1(i,zz)-v3d(i,k,j))*rdelt\r
327           pek(i,k,j)     = pek2d(i,zz)\r
328           pep(i,k,j)     = pep2d(i,zz)\r
329         enddo\r
330         pp = pp + 1\r
331       enddo\r
333       pp = 0\r
334       do k=kts,kte+1\r
335         zz = kte+1-pp\r
336         do i=its,ite\r
337           exch_h(i,k,j)  = exh2d(i,zz)\r
338           exch_m(i,k,j)  = exm2d(i,zz)\r
339         enddo\r
340         pp = pp + 1\r
341       enddo\r
343       do i=its,ite\r
344         pblh(i,j) = pblh1(i)\r
345         kpbl(i,j) = kx-kpbl1(i)+1\r
346       end do\r
348       if(bl_eeps_tkeadv==1)then\r
349         do k=kts,kte\r
350           do i=its,ite\r
351              pek_adv(i,k,j) = pek(i,k,j)\r
352              pep_adv(i,k,j) = pep(i,k,j)\r
353           end do\r
354         end do\r
355       end if\r
357    enddo\r
360    end subroutine eeps\r
362 !-------------------------------------------------------------------------------\r
364    subroutine eeps2d(pu,pv,tz,pqv,pqc,pqi,pqr,pqs,pqg,prs,poz,zz,rho,rthraten,pek,pep,       &\r
365                   psfcpa,ust,rmol,wspd,xland,qsfc,hfx,qfx,tsk,                  &\r
366                   dt,dx,itimestep,k_h,k_m,pblh,kpbl,lq,km)\r
367 !-------------------------------------------------------------------------------\r
368    implicit none\r
369 !-------------------------------------------------------------------------------\r
371    integer,  intent(in   )   ::     itimestep,lq,km\r
372    real,     intent(in   )   ::     dt, dx\r
374    real,     dimension( lq )                                            , &\r
375              intent(in)   ::                                              ust, &\r
376                                                                         xland, &\r
377                                                                           hfx, &\r
378                                                                           qfx, &\r
379                                                                           rmol,&\r
380                                                                           wspd,&\r
381                                                                    psfcpa,qsfc,&\r
382                                                                            tsk\r
384    real,     dimension( lq,km ),intent(inout)      ::          pu, &\r
385                                                                pv, &\r
386                                                                tz, &\r
387                                                               pqv, &\r
388                                                               pqc, &\r
389                                                               pqi, &\r
390                                                               pek, &\r
391                                                               pep, &\r
392                                                               rho, &\r
393                                                          rthraten\r
394   \r
395    real,     dimension( lq,km ),intent(in)      ::            prs, &\r
396                                                               poz, &\r
397                                                               pqr, &\r
398                                                               pqs, &\r
399                                                               pqg\r
400    real,     dimension( lq,km+1 ),intent(in)               ::  zz\r
402    real,     dimension( lq,km+1 ),intent(out)              :: k_h,k_m\r
403    real,     dimension( lq ),intent(out)                   ::  pblh    \r
404    integer,  dimension( lq ),intent(out)                   ::  kpbl\r
406 ! Local Vars\r
407    real,     dimension( lq,km ) :: fac,psp,pt,ztv,szv,disht,up2,cpm\r
408    real,     dimension( lq,km ) :: qvs,st,sh,rich,richb,dum2,dzz,doz,st0,rich0\r
409    real,     dimension( lq,km ) :: am,src,src1,ax,bx,cx,ac,ab,ba,alt\r
410    real,     dimension( lq,km-1 ) :: ax1,bx1,cx1,yy1,amt\r
411    real,     dimension( lq,km+1 ) :: akm,akh,pnt\r
413    real,     dimension( lq)    :: tkm,zkm,qkm,us2,us,dens,wstar2,sflux\r
414    real,     dimension( lq)    :: rm2,sfa,hs,hq\r
415    real,     dimension( lq)    :: gamat,gamaq\r
416    real,     dimension( lq,km) :: sgk1,sgk2,edrain \r
417    integer ::  j,k,mdt,lvl\r
418    real    ::  dtt,dum,dum0,alh,eqp,qvsw,qvsi,faf,rdz,du,dv,dss,bvf,rbdts\r
419    real    ::  tk,qvk,alk,af,ag,aks,aco\r
420    real    ::  duj,alv0,richf,rch,alv,ak\r
421    real    ::  govrth,bfx0,wstar3,coef\r
423 ! Local Vars for top-down diffusion\r
424    integer,  dimension( lq)    :: kminrad,kcld,k700,k950\r
425    logical,  dimension( lq)    :: flg,flg1,scuflg\r
426    real,     dimension( lq)    :: zl1,wm3,ent_eff,radsum,minrad,zminrad\r
427    real,     dimension( lq)    :: eis\r
428    real,     dimension( lq,km) :: pci,TKEprodTD,zfac,zfacent\r
429    real    :: dthvx,tmp1,temps,templ\r
430    integer :: ktop,kk1,kk2\r
431    real    :: radflux,rcldb,rvls,he0,he1,he2\r
432 !end top down\r
433    \r
434    mdt  = 1 + int(sqrt(dt/9.0)+0.00001)\r
435    dtt = dt/mdt\r
437 !  dzz for dz between full levels, doz for dz between half levels\r
438    do k=1,km\r
439      do j=1,lq\r
440        dzz(j,k) = zz(j,k)-zz(j,k+1)\r
441        alt(j,k) = dt/dzz(j,k)\r
442        up2(j,k)=grv/max(us2min,pu(j,k)*pu(j,k)+pv(j,k)*pv(j,k))\r
443        TKEprodTD(j,k) = 0.0\r
444        zfac(j,k) = 0.0\r
445        zfacent(j,k) = 0.0\r
446      enddo\r
447    enddo\r
449    do k=1,km-1\r
450      do j=1,lq\r
451        doz(j,k)  = poz(j,k)-poz(j,k+1)\r
452      end do\r
453    end do\r
455    do j=1,lq\r
456      doz(j,km)= 2.0*poz(j,km) ! not used\r
457    enddo\r
459    do k=1,km-1\r
460      do j=1,lq\r
461        amt(j,k) = dtt/doz(j,k)\r
462      enddo\r
463    enddo\r
465    do k=1,km\r
466      do j=1,lq\r
467        sgk1(j,k) = 0.1\r
468        sgk2(j,k) = 0.01\r
469      end do\r
470    end do\r
472    do k=1,km\r
473    do j=1,lq\r
474       psp(j,k)=(prs(j,k)/p1000mb)**rcp\r
475       pt(j,k)=tz(j,k)/psp(j,k)\r
476       cpm(j,k) = cp*(1.+0.8*pqv(j,k))\r
477       dum=1.0+ep1*pqv(j,k)\r
478       ztv(j,k)=tz(j,k)*dum\r
479       szv(j,k)=pt(j,k)*dum\r
480       pci(j,k) = pqc(j,k) + pqi(j,k)\r
481       fac(j,k) = pqc(j,k) + pqi(j,k) + pqr(j,k) + pqs(j,k) + pqg(j,k)\r
482       cx(j,k)=1.0+pqv(j,k)+fac(j,k)\r
484       alh=3.1484e6-2.37e3*tz(j,k)\r
485       eqp=611.2*exp(cs2*(tz(j,k)-cs3)/(tz(j,k)-cs4))\r
486       eqp=min(0.5*prs(j,k),eqp)\r
487       qvsw=0.622*eqp/(prs(j,k)-eqp)\r
488       eqp=611.0*exp(ci2*(tz(j,k)-ci3)/(tz(j,k)-ci4))\r
489       eqp=min(0.5*prs(j,k),eqp)\r
490       qvsi=0.622*eqp/(prs(j,k)-eqp)\r
491       if(tz(j,k).ge.t0) then\r
492         faf=0.0\r
493       else if(tz(j,k).le.hgfr) then\r
494         faf=1.0\r
495       else\r
496         if(pqi(j,k).le.epsl) then\r
497         faf=0.0\r
498         else if(pqc(j,k).le.epsl) then\r
499         faf=1.0\r
500         else\r
501         faf=pqi(j,k)/(pqc(j,k)+pqi(j,k))\r
502         endif\r
503       endif\r
504       qvs(j,k)=(1.0-faf)*qvsw+faf*qvsi\r
505       ax(j,k)=(1.0-faf)*alh+faf*xls\r
506    end do\r
507    end do\r
510 !  to calculate the Richardson number and vertical shear of\r
511 !  the horizontal wind above the surface layer\r
513       do k=1,km-1\r
514       do j=1,lq\r
515         rdz=1.0/(poz(j,k)-poz(j,k+1))\r
516         du=pu(j,k)-pu(j,k+1)\r
517         dv=pv(j,k)-pv(j,k+1)\r
518         sh(j,k)=(du*du+dv*dv)*(rdz*rdz)\r
519         if((pqv(j,k).lt.qvs(j,k)).or.(pqv(j,k+1).lt.qvs(j,k+1))) then\r
520           dss=log(szv(j,k)/szv(j,k+1))\r
521           st(j,k)=grv*dss*rdz\r
522         else\r
523           dss=log(pt(j,k)/pt(j,k+1))\r
524           tk=0.5*(tz(j,k+1)+tz(j,k))\r
525           qvk=0.5*(pqv(j,k+1)+pqv(j,k))\r
526           alk=0.5*(ax(j,k+1)+ax(j,k))\r
527           af=alk*qvk/(rgas*tk)\r
528           ag=alk/(cp*tk)\r
529           aks=ep2*af*ag\r
530           aco=(1.0+af)/(1.0+aks)\r
531           st(j,k)=(aco*(dss+ag*(pqv(j,k)-pqv(j,k+1))) &\r
532                -(cx(j,k)-cx(j,k+1)))*grv*rdz\r
533         endif\r
534         rich(j,k)=st(j,k)/max(1.0e-7,sh(j,k))\r
535       end do\r
536       end do\r
538 ! calculate first guess pblh\r
539       do k=1,km\r
540       do j=1,lq\r
541         richb(j,k)=poz(j,k)*(szv(j,k)/szv(j,km)-1.0)*up2(j,k)\r
542       enddo\r
543       enddo\r
545       do j=1,lq\r
546         pblh(j)=poz(j,km)\r
547         kpbl(j)=km\r
548         lvl=km\r
549         do k=km-1,5,-1\r
550           if(richb(j,k).ge.ricr) then\r
551             lvl=k\r
552             exit\r
553           endif\r
554         enddo\r
556         if(lvl.lt.km) then\r
557           dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
558           pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
559           pblh(j)=min(pblh(j),poz(j,lvl))\r
560           kpbl(j)=lvl\r
561         endif\r
562       end do\r
564 !  preparing for tke/ep calculations\r
566       do j=1,lq\r
567         tkm(j) = tz(j,km)\r
568         zkm(j) = poz(j,km)\r
569         qkm(j) = pqv(j,km)\r
570         sfa(j) = tsk(j)*(p1000mb/psfcpa(j))**rcp\r
571         dens(j) = psfcpa(j)/(rgas*tkm(j)*(1.+ep1*qkm(j)))\r
573         sh(j,km)=0.0\r
574         st(j,km)=0.0\r
575         rm2(j)=ust(j)*ust(j)/wspd(j)\r
576         hs(j)=hfx(j)/(dens(j)*cpm(j,km))\r
577         hq(j)=qfx(j)/dens(j)\r
579         duj=ust(j)*ust(j)\r
580         pep(j,km)=duj*ust(j)/(akv*zkm(j))\r
581         pep(j,km) = min(maxpep,max(minpep,pep(j,km)))\r
582         govrth = grv/pt(j,km)\r
583         sflux(j) = hs(j) + hq(j)*ep1*pt(j,km)\r
584         bfx0    = max(sflux(j),0.)\r
585         wstar3  = govrth*bfx0*pblh(j)\r
586         wstar2(j)  = (wstar3**2)**t13\r
588         rich(j,km) = grv*zkm(j)*  &\r
589             log(szv(j,km)/(sfa(j)*(1.+ep1*qsfc(j))))/(wspd(j)**2.)\r
591         if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
592           pek(j,km)=3.75*duj+0.2*wstar2(j)+duj*((-zkm(j)*rmol(j))**2.)**t13\r
593         else\r
594           pek(j,km)=3.75*duj\r
595         endif\r
596         pek(j,km) = min(maxpek,max(minpek,pek(j,km)))\r
598       end do\r
600 ! to recalculate pblh here\r
601       do j=1,lq\r
602         if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
603           dss=szv(j,km)+min(5.0,bt*sflux(j)/(max(0.1,sqrt(wstar2(j)))))\r
604           do k=1,km\r
605             richb(j,k)=poz(j,k)*(szv(j,k)/dss-1.0)*up2(j,k)\r
606           end do\r
608           lvl=km\r
609           do k=km-1,5,-1\r
610             if(richb(j,k).ge.ricr) then\r
611               lvl=k\r
612               exit\r
613             end if\r
614           enddo\r
615  \r
616           if(lvl.lt.km) then\r
617             dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
618             pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
619             pblh(j)=min(pblh(j),poz(j,lvl))\r
620             kpbl(j)=lvl\r
621           end if\r
622         end if\r
623       end do\r
625 ! to add the countergradient term\r
626       do j=1,lq\r
627         if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
628           du=max(0.1,sqrt(wstar2(j)))\r
629           gamat(j)=max(bt*hs(j)/(pblh(j)*du),0.)\r
630           gamaq(j)=max(bt*hq(j)/(pblh(j)*du),0.)\r
631         else\r
632           gamat(j)=0.0\r
633           gamaq(j)=0.0\r
634         endif\r
635       end do\r
637 ! to add TKE source driven by stratocumulus cloud top cooling\r
638 ! a switch bl_eeps_topdown controls if it is on or off\r
639     if (bl_eeps_topdown.eq.1)then\r
640       do j=1,lq\r
641          scuflg(j)=.true.\r
642          flg(j)   =.true.\r
643          flg1(j)  =.true.\r
644          minrad(j)  = 100.\r
645          kminrad(j) = km\r
646          kcld(j)    = km\r
647          k700(j)    = km\r
648          k950(j)    = km\r
649          zminrad(j) = poz(j,km)\r
650          zl1(j)     = poz(j,km)\r
651       end do\r
653 !  look for stratocumulus\r
654 !  look for the highest model level we will use\r
655       do k=km,2,-1\r
656       do j=1,lq\r
657         if(flg(j) .and. (prs(j,km)-prs(j,k))>3.0e4) then\r
658           k700(j)=k\r
659           flg(j)=.false.\r
660         end if\r
661         if(flg1(j) .and. (prs(j,km)-prs(j,k))>5.e3) then\r
662           k950(j)=k\r
663           flg1(j)=.false.\r
664         end if\r
665       end do\r
666       end do\r
667       ktop = minval(k700)\r
669 !  calculate EIS based on Marquet and Bechtold (2020)\r
670       do j=1,lq\r
671           kk1 = k700(j)\r
672           kk2 = k950(j)\r
673           he1 = cp*tz(j,kk1)*(1.+5.87*(pqv(j,kk1)+pci(j,kk1))) &\r
674                -xlv*pqc(j,kk1)-xls*pqi(j,kk1)+grv*poz(j,kk1)\r
675           he2 = cp*tz(j,kk2)*(1.+5.87*(pqv(j,kk2)+pci(j,kk2))) &\r
676                -xlv*pqc(j,kk2)-xls*pqi(j,kk2)+grv*poz(j,kk2)\r
677           he0 = cp*tz(j,km)*(1.+5.87*(pqv(j,km)+pci(j,km))) &\r
678                -xlv*pqc(j,km)-xls*pqi(j,km)+grv*poz(j,km) \r
679           eis(j) = max(he1-he2,he2-he0)\r
680           eis(j) = eis(j)/cp\r
681           if(eis(j) < 5) scuflg(j)=.false.\r
682       end do\r
684 !  look for the highest cloud level\r
685       do j=1,lq\r
686          flg(j)=scuflg(j)\r
687       end do\r
689       do k=ktop,km\r
690       do j=1,lq\r
691         if(flg(j) .and. k >= k700(j) .and. pci(j,k) >= 1.e-6) then\r
692            kcld(j)=k\r
693            flg(j)=.false.\r
694         end if\r
695       end do\r
696       end do\r
698 !   find the level of the minimum radiative cooling\r
699       do k=ktop,km\r
700       do j=1,lq\r
701         if(scuflg(j)) then\r
702           if(k >= kcld(j) .and. pci(j,k) >= 1.e-6) then\r
703             if(rthraten(j,k) < minrad(j)) then\r
704               minrad(j)=rthraten(j,k)\r
705               kminrad(j)=k\r
706               zminrad(j)=poz(j,k)\r
707             endif\r
708           endif\r
709         endif\r
710       end do\r
711       end do\r
713 ! to exclude some scenarios\r
714 !  1) cloud top should be below PBL height\r
715 !  2) minimum radiative cooling level can't too close to surface\r
716 !  3) minimum radiative cooling should be negative\r
717       do j=1,lq\r
718         if(scuflg(j) .and. minrad(j)>=0.)              scuflg(j)=.false.\r
719         if(scuflg(j) .and. kcld(j)<kpbl(j))            scuflg(j)=.false.\r
720         if(scuflg(j) .and. kminrad(j)>=(km-1))         scuflg(j)=.false.\r
721       end do\r
723       do j=1,lq\r
724              if (scuflg(j)) then\r
725                 k = kminrad(j)\r
726                 templ=tz(j,k)\r
727                 !rvls is ws at full level\r
728                 rvls=100.*6.112*exp(17.67*(templ-273.16)/(templ-29.65))*(ep2/prs(j,k))\r
729                 temps=templ + (pqv(j,k)+pci(j,k)-rvls)/(cp/xlv  +  ep2*xlv*rvls/(rgas*templ**2))\r
730                 rvls=100.*6.112*exp(17.67*(temps-273.15)/(temps-29.65))*(ep2/prs(j,k))\r
731                 rcldb=max(pqv(j,k)+pci(j,k)-rvls,0.)\r
733                 !entrainment efficiency\r
734                 dthvx     = szv(j,k-2)-szv(j,k)\r
735                 dthvx     = max(dthvx,0.1)\r
736                 tmp1      = xlv / cp * rcldb/(psp(j,k)*dthvx)\r
737                 !Originally from Nichols and Turton (1986), where a2 = 60, but lowered\r
738                 !here to 8, as in Grenier and Bretherton (2001).\r
739                 ent_eff(j)   = 0.2 + 0.2*8.*tmp1\r
740              end if\r
741       end do\r
743       do k=km,ktop,-1\r
744          do j=1,lq\r
745              if(scuflg(j) .and. k<=(kminrad(j)+1) .and. k>=kcld(j)) then\r
746                    radflux=rthraten(j,k)*psp(j,k)         !converts theta/s to temp/s\r
747                    radflux=radflux*cp/grv*(prs(j,k)-prs(j,k-1)) ! converts temp/s to W/m^2\r
748                    if (radflux < 0.0 ) radsum(j)=abs(radflux)+radsum(j)\r
749              end if\r
750          end do\r
751       end do\r
753       do j=1,lq\r
754              if(scuflg(j)) then\r
755                 k = kminrad(j)\r
756                 radsum(j)=min(radsum(j),120.0)\r
757                 bfx0 = max(radsum(j)/rho(j,k)/cp,0.)\r
758                 wm3(j)    = grv/szv(j,k)*bfx0*min(pblh(j),1500.) ! this is wstar3(i)\r
759              end if\r
760       end do\r
762       do k=km,ktop,-1\r
763          do j=1,lq\r
764              if(scuflg(j) .and. k>=kcld(j)) then\r
765                    !Analytic vertical profile\r
766                    zfac(j,k) = min(max((1.-(poz(j,k)-zl1(j))/(zminrad(j)-zl1(j))),zfmin),1.)\r
767                    zfacent(j,k) = 10.*max((zminrad(j)-poz(j,k))/zminrad(j),0.0)*(1.-zfac(j,k))**3\r
768                    !Calculate TKE production = 2(g/TH)(w'TH'), where w'TH' = A(TH/g)wstar^3/PBLH,\r
769                    !A = ent_eff, and wstar is associated with the radiative cooling at top of PBL.\r
770                    !An analytic profile controls the magnitude of this TKE prod in the vertical. \r
771                    TKEprodTD(j,k)=2.*ent_eff(j)*wm3(j)/max(pblh(j),100.)*zfacent(j,k)\r
772                    TKEprodTD(j,k)= max(TKEprodTD(j,k),0.0)\r
773              end if !end cloud check\r
774           end do ! end lq loop\r
775        end do  ! end km loop\r
777     end if !end top-down check\r
779 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
780 !  time integration of tke and ep equations\r
781 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
782       do k=1,km\r
783       do j=1,lq\r
784         richf=0.1620092\r
785         if(rich(j,k).lt.0.145) then\r
786           richf=0.6588*(rich(j,k)+0.1776-sqrt(rich(j,k)*rich(j,k)- &\r
787               0.3221*rich(j,k)+0.03156))\r
788           richf=max(-1000.0,min(0.1620092,richf))\r
789         endif\r
790         dum2(j,k)=1.116893\r
791         if(richf.lt.0.1620092) dum2(j,k)=1.318*(0.2231-richf) &\r
792               /(0.2341-richf)\r
793       end do\r
794       end do\r
796       lvl=1\r
797       do while(lvl .le. mdt)\r
799       do k=1,km-1\r
800         do j=1,lq\r
801           akm(j,k)=c2*pek(j,k)*pek(j,k)/pep(j,k)\r
802           akh(j,k)=dum2(j,k)*akm(j,k)\r
803           akm(j,k)=sgk1(j,k)+akm(j,k)\r
804           akh(j,k)=sgk2(j,k)+akh(j,k)\r
805           akm(j,k)=min(akvmx,akm(j,k))\r
806           akh(j,k)=min(akvmx,akh(j,k))\r
807           if(st(j,k).gt.0.0) then\r
808             bvf=sqrt(st(j,k))\r
809             rbdts=min(1.,sqrt(rich(j,k)/0.145))\r
810             edrain(j,k)=0.1/c3*rbdts*bvf*pek(j,k)\r
811           else\r
812             edrain(j,k)=0.0 \r
813           endif\r
814           if(poz(j,k).gt.pblh(j))then\r
815             src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*st(j,k)+TKEprodTD(j,k)\r
816             src1(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*st(j,k)+edrain(j,k)+TKEprodTD(j,k)\r
817           else\r
818             src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
819                    -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k)))+TKEprodTD(j,k)\r
820             src1(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
821                    -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k)))+edrain(j,k)+TKEprodTD(j,k)\r
822           end if\r
823         end do\r
824       end do\r
826       do j=1,lq\r
827         akm(j,km)=c2*pek(j,km)*pek(j,km)/pep(j,km)\r
828         akh(j,km)=dum2(j,km)*akm(j,km)\r
829       end do\r
831 !  to calculate the vertical diffusion coefficients at u,v,levels\r
833       do k=2,km\r
834       do j=1,lq\r
835         am(j,k)=0.5*(akm(j,k-1)+akm(j,k))/dzz(j,k)\r
836       end do\r
837       end do\r
839       do j=1,lq\r
840         am(j,1)=0.5*akm(j,1)/dzz(j,1)\r
841       end do\r
843 !  to calculate the cofficients for the triangle matrix\r
844 !  for tke dissipation and solve the matrix\r
846       do k=1,km-1\r
847         do j=1,lq\r
848           dum=(c3*src1(j,k)-c4*pep(j,k))*dtt/pek(j,k)\r
849           ax1(j,k)=-c5*amt(j,k)*am(j,k+1)\r
850           cx1(j,k)=-c5*amt(j,k)*am(j,k)\r
851           if(dum.le.0.5) then\r
852             bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))-dum\r
853             yy1(j,k)=pep(j,k)\r
854           else\r
855             bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))\r
856             yy1(j,k)=pep(j,k)+dum*pep(j,k)\r
857           endif\r
858         end do\r
859       end do\r
861       do j=1,lq\r
862         cx1(j,1)=0.0\r
863         ax1(j,km-1)=0.0\r
864         yy1(j,km-1)=yy1(j,km-1)+c5*amt(j,km-1)*am(j,km)*pep(j,km)\r
865       end do\r
867       call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
869 !  to calculate the cofficients for the triangle matrix\r
870 !  for tke and solve the matrix\r
872       do k=1,km-1\r
873         do j=1,lq\r
874           pep(j,k)=min(maxpep,max(minpep, yy1(j,k)))\r
875           dum=(src(j,k)-pep(j,k))*dtt\r
876           ax1(j,k)=-c1*amt(j,k)*am(j,k+1)\r
877           bx1(j,k)=1.0+c1*amt(j,k)*(am(j,k)+am(j,k+1))\r
878           cx1(j,k)=-c1*amt(j,k)*am(j,k)\r
879           yy1(j,k)=pek(j,k)+dum\r
880         end do\r
881       end do\r
883       do j=1,lq\r
884         cx1(j,1)=0.0\r
885         ax1(j,km-1)=0.0\r
886         yy1(j,km-1)=yy1(j,km-1)+c1*amt(j,km-1)*am(j,km)*pek(j,km)\r
887       end do\r
889       call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
891       do k=1,km-1\r
892       do j=1,lq\r
893         pek(j,k)=min(maxpek,max(minpek, yy1(j,k)))\r
894       end do\r
895       end do\r
897       lvl = lvl + 1\r
898       end do ! end while\r
900 !-------------------------------------------------------------\r
901 !  to calculate the vertical diffusion coeficient for momentum \r
902 !  heat, moisture\r
903 !-------------------------------------------------------------\r
904 !  note: akm and akh at the first level (km) will not be used \r
905       do k=2,km\r
906         do j=1,lq\r
907           akm(j,k)=c2*pek(j,k-1)*pek(j,k-1)/pep(j,k-1)\r
908           akh(j,k)=dum2(j,k)*akm(j,k)\r
909           akm(j,k)=sgk1(j,k)+akm(j,k)\r
910           akh(j,k)=sgk2(j,k)+akh(j,k)\r
911           akm(j,k)=min(akvmx,akm(j,k))\r
912           akh(j,k)=min(akvmx,akh(j,k))\r
913           pnt(j,k)=akh(j,k)\r
914           if(poz(j,k).gt.pblh(j)) pnt(j,k)=0.\r
915           akm(j,k)=akm(j,k)/doz(j,k-1)\r
916           akh(j,k)=akh(j,k)/doz(j,k-1)\r
917         end do\r
918       end do\r
920       do j=1,lq\r
921         akm(j,1)=0.0\r
922         akh(j,1)=0.0\r
923         pnt(j,1)=0.0\r
924         akm(j,km+1)=0.0\r
925         akh(j,km+1)=0.0\r
926         pnt(j,km+1)=0.0\r
927       end do\r
929 ! Momentum fluxes calculation\r
931       do k=1,km\r
932       do j=1,lq\r
933         ax(j,k)=-alt(j,k)*akm(j,k+1)\r
934         bx(j,k)=1.0+alt(j,k)*(akm(j,k)+akm(j,k+1))\r
935         cx(j,k)=-alt(j,k)*akm(j,k)\r
936         ba(j,k)=ax(j,k)\r
937         ab(j,k)=bx(j,k)\r
938         ac(j,k)=cx(j,k)\r
939       end do\r
940       end do\r
942       do j=1,lq\r
943         bx(j,km) = bx(j,km)+alt(j,km)*rm2(j)\r
944         ab(j,km)=bx(j,km)\r
945       end do\r
948       call tridiag(ba,ab,ac,pu,lq,km)\r
950       call tridiag(ax,bx,cx,pv,lq,km)\r
952 ! Heat and moisture fluxes calculation\r
954       do k=1,km\r
955       do j=1,lq\r
956         ax(j,k)=-alt(j,k)*akh(j,k+1)\r
957         bx(j,k)=1.0+alt(j,k)*(akh(j,k)+akh(j,k+1))\r
958         cx(j,k)=-alt(j,k)*akh(j,k)\r
959         ba(j,k)=ax(j,k)\r
960         ab(j,k)=bx(j,k)\r
961         ac(j,k)=cx(j,k)\r
962         dum0=0.0       \r
963         if(k.gt.1) dum0=pep(j,k-1)\r
964         dum=pep(j,k)\r
965         if(dum0.le.1.e-6) dum0=0.0\r
966         if(dum.le.1.e-6) dum=0.0\r
967         disht(j,k)=0.5*(dum0+dum)/cpm(j,k)\r
968         pt(j,k)=pt(j,k)+disht(j,k)*dt/psp(j,k)+alt(j,k)* &\r
969                 gamat(j)*(pnt(j,k+1)-pnt(j,k))\r
970       end do\r
971       end do\r
973       do j=1,lq\r
974         pt(j,km)=pt(j,km)+alt(j,km)*hs(j)\r
975       end do\r
977       call tridiag(ba,ab,ac,pt,lq,km)\r
979       do k=1,km\r
980       do j=1,lq\r
981         ba(j,k)=ax(j,k)\r
982         ab(j,k)=bx(j,k)\r
983         ac(j,k)=cx(j,k)\r
984       end do\r
985       end do\r
987       do k=1,km\r
988       do j=1,lq\r
989          pqv(j,k)=pqv(j,k)+alt(j,k)*gamaq(j)* &\r
990                  (pnt(j,k+1)-pnt(j,k))\r
991       end do\r
992       end do\r
994       do j=1,lq\r
995         pqv(j,km)=pqv(j,km)+alt(j,km)*hq(j)\r
996       end do\r
998       call tridiag(ba,ab,ac,pqv,lq,km)\r
1000       do k=1,km\r
1001       do j=1,lq\r
1002         ba(j,k)=ax(j,k)\r
1003         ab(j,k)=bx(j,k)\r
1004         ac(j,k)=cx(j,k)\r
1005       end do\r
1006       end do\r
1008       call tridiag(ba,ab,ac,pqc,lq,km)\r
1010       call tridiag(ax,bx,cx,pqi,lq,km)\r
1012       do k=1,km\r
1013       do j=1,lq\r
1014         tz(j,k)=pt(j,k)*psp(j,k)\r
1015         pqv(j,k)=max(epsl,pqv(j,k))\r
1016         pqc(j,k)=max(epsl,pqc(j,k))\r
1017         pqi(j,k)=max(epsl,pqi(j,k))\r
1018 !-----------------------------------------------------\r
1019 !  immediate melting of cloud ice below freezing level\r
1020 !-----------------------------------------------------\r
1021         if(tz(j,k).ge.t0.and.pqi(j,k).gt.epsl) then\r
1022           if(pqc(j,k).le.epsl) then\r
1023             pqc(j,k)=pqi(j,k)\r
1024           else\r
1025             pqc(j,k)=pqc(j,k)+pqi(j,k)\r
1026           endif\r
1027           alh=3.1484e6-2.37e3*tz(j,k)\r
1028           tz(j,k)=tz(j,k)-(xls-alh)*pqi(j,k)/cpm(j,k)\r
1029           pqi(j,k)=epsl\r
1030         endif\r
1031       end do\r
1032       end do\r
1034       do k=2,km\r
1035       do j=1,lq\r
1036         k_h(j,k)=akh(j,k)*doz(j,k-1)\r
1037         k_m(j,k)=akm(j,k)*doz(j,k-1)\r
1038       end do\r
1039       end do\r
1041       do j=1,lq\r
1042         k_h(j,1)=0.\r
1043         k_m(j,1)=0.\r
1044         k_m(j,km+1)= 0. !c2*pek(j,km)*pek(j,km)/pep(j,km) ! we don't need it\r
1045         k_h(j,km+1)= 0. !dum2(j,km)*k_m(j,km+1)           ! we don't need it\r
1046       end do\r
1048    end subroutine eeps2d\r
1050 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
1052 ! ==================================================================\r
1053    subroutine tridiag(a,b,c,d,lq,km)\r
1054 !! to solve system of linear eqs on tridiagonal matrix n times n\r
1055 !! after Peaceman and Rachford, 1955\r
1056 !! a,b,c,d - are vectors of order n \r
1057 !! a,b,c - are coefficients on the LHS\r
1058 !! d - is initially RHS on the output becomes a solution vector\r
1060 !-------------------------------------------------------------------\r
1061     implicit none\r
1062     INTEGER, INTENT(in):: lq,km\r
1063     REAL, DIMENSION(lq,km), INTENT(in) :: a,b\r
1064     REAL, DIMENSION(lq,km), INTENT(inout) :: c,d\r
1066     INTEGER :: j,k\r
1067     REAL :: p\r
1068     REAL, DIMENSION(lq,km) :: q\r
1070     do j=1,lq\r
1071       c(j,1)=0.\r
1072       q(j,km)=-c(j,km)/b(j,km)\r
1073       d(j,km)=d(j,km)/b(j,km)\r
1074     end do\r
1076     DO k=km-1,1,-1\r
1077     do j=1,lq\r
1078        p=1./(b(j,k)+a(j,k)*q(j,k+1))\r
1079        q(j,k)=-c(j,k)*p\r
1080        d(j,k)=(d(j,k)-a(j,k)*d(j,k+1))*p\r
1081     end do\r
1082     ENDDO\r
1084     DO k=2,km\r
1085     do j=1,lq\r
1086        d(j,k)=d(j,k)+q(j,k)*d(j,k-1)\r
1087     end do\r
1088     ENDDO\r
1090    end subroutine tridiag\r
1092 !-------------------------------------------------------------------------------\r
1093    subroutine eepsinit(rublten,rvblten,rthblten,rqvblten,                       &\r
1094                       rqcblten,rqiblten,p_qi,p_first_scalar,pek,pep,            &\r
1095                       restart, allowed_to_read,                                &\r
1096                       ids, ide, jds, jde, kds, kde,                            &\r
1097                       ims, ime, jms, jme, kms, kme,                            &\r
1098                       its, ite, jts, jte, kts, kte                 )\r
1099 !-------------------------------------------------------------------------------\r
1100    implicit none\r
1101 !-------------------------------------------------------------------------------\r
1103    logical , intent(in)          :: restart, allowed_to_read\r
1104    integer , intent(in)          ::  ids, ide, jds, jde, kds, kde,             &\r
1105                                      ims, ime, jms, jme, kms, kme,             &\r
1106                                      its, ite, jts, jte, kts, kte\r
1107    integer , intent(in)          ::  p_qi,p_first_scalar\r
1108    real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) ::             &\r
1109                                                                       rublten, &\r
1110                                                                       rvblten, &\r
1111                                                                      rthblten, &\r
1112                                                                      rqvblten, &\r
1113                                                                      rqcblten, &\r
1114                                                                      rqiblten\r
1115    real, dimension( ims:ime, kms:kme, jms:jme ),intent(out)   ::     pek, pep\r
1116    integer :: i, j, k, itf, jtf, ktf\r
1118    jtf = min0(jte,jde-1)\r
1119    ktf = min0(kte,kde-1)\r
1120    itf = min0(ite,ide-1)\r
1122    if(.not.restart)then\r
1123      do j = jts,jtf\r
1124        do k = kts,ktf\r
1125          do i = its,itf\r
1126             rublten(i,k,j) = 0.\r
1127             rvblten(i,k,j) = 0.\r
1128             rthblten(i,k,j) = 0.\r
1129             rqvblten(i,k,j) = 0.\r
1130             rqcblten(i,k,j) = 0.\r
1131             pek(i,k,j) = minpek\r
1132             pep(i,k,j) = minpep\r
1133          enddo\r
1134        enddo\r
1135      enddo\r
1136    endif\r
1138    if (p_qi .ge. p_first_scalar .and. .not.restart) then\r
1139      do j = jts,jtf\r
1140        do k = kts,ktf\r
1141          do i = its,itf\r
1142            rqiblten(i,k,j) = 0.\r
1143          enddo\r
1144        enddo\r
1145      enddo\r
1146    endif\r
1148    end subroutine eepsinit\r
1149 !-------------------------------------------------------------------------------  \r
1150 end module module_bl_eeps\r