37 subroutine sgscloud_radpre_run( &
39 flag_init,flag_restart, &
40 con_g, con_pi, eps, epsm1, &
44 qc, qi, qv, T3D, P3D, exner, &
46 qci_conv,qlc,qli,ud_mf, &
47 imfdeepcnv, imfdeepcnv_gf, &
50 qc_save, qi_save, qs_save, &
51 qc_bl,qi_bl,cldfra_bl, &
52 delp,clouds1,clouds2,clouds3, &
54 clouds8,clouds9,slmsk, &
55 nlay, plyr, xlat, dz,de_lgth, &
57 imp_physics, imp_physics_gfdl,&
70 real(kind=kind_phys),
intent(in) :: con_g, con_pi, eps, epsm1
71 real(kind=kind_phys),
intent(in) :: r_v, cpv, rcp
72 real(kind=kind_phys),
intent(in) :: xlv, xlf, cp
73 real(kind=kind_phys),
intent(in) :: dt,fhswr
74 real :: xls, xlvcp, xlscp
75 real(kind=kind_phys) ::
gfac
76 integer,
intent(in) :: im, levs, imfdeepcnv, imfdeepcnv_gf, &
77 & nlay, imfdeepcnv_sas, imfdeepcnv_c3, imp_physics, &
78 & imp_physics_gfdl, imp_physics_fa
79 logical,
intent(in) :: flag_init, flag_restart, do_mynnedmf
81 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: qc, qi
82 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: qr, qs, qg
84 real(kind=kind_phys),
dimension(:,:),
intent(inout),
optional :: qci_conv
85 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: qlc, qli
86 real(kind=kind_phys),
dimension(:,:),
intent(in),
optional :: ud_mf
87 real(kind=kind_phys),
dimension(:,:),
intent(in) :: t3d,delp
88 real(kind=kind_phys),
dimension(:,:),
intent(in) :: qv,p3d,exner
89 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: &
90 & clouds1,clouds2,clouds3,clouds4,clouds5, &
92 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: qc_save, qi_save, qs_save
93 real(kind=kind_phys),
dimension(:,:),
intent(in),
optional :: qc_bl, qi_bl, cldfra_bl
94 real(kind=kind_phys),
dimension(:),
intent(in) :: slmsk, xlat, de_lgth
95 real(kind=kind_phys),
dimension(:,:),
intent(in) :: plyr, dz
96 real(kind=kind_phys),
dimension(:,:),
intent(inout) :: cldsa
97 integer,
dimension(:,:),
intent(inout) :: mbota, mtopa
98 integer,
intent(in) :: iovr
99 character(len=*),
intent(out) :: errmsg
100 integer,
intent(out) :: errflg
103 real(kind=kind_phys) :: ptop1(im,3+1)
104 real(kind=kind_phys) ::
ptopc(3+1,2 )
106 data ptopc / 1050., 650., 400., 0.0, 1050., 750., 500., 0.0 /
107 real(kind=kind_phys),
dimension(im,nlay) :: cldcnv
108 real(kind=kind_phys),
dimension(im) :: rxlat
109 real(kind=kind_phys) :: tc, tk, liqfrac, iwc, ice_frac, snow_frac
112 real(kind=kind_phys),
dimension(im,nlay) :: alpha_dummy
116 real,
parameter :: coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100.
117 real :: rhgrid,h2oliq,qsat,tem1,tem2,clwt,es,onemrh,value
120 real :: a, f, sigq, qmq, qt, xl, th, thl, rsl, cpm, cb_cf
121 real(kind=kind_phys) :: tlk
124 integer,
parameter :: conv_cf_opt = 0
138 if (flag_init .and. (.not. flag_restart))
then
147 if ( qi(i,k) > 1e-7 .OR. qc(i,k) > 1e-7 )
then
148 es = min( p3d(i,k), fpvs( t3d(i,k) ) )
149 qsat = max( qmin, eps * es / (p3d(i,k) + epsm1*es) )
150 rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
151 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k)
152 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
154 if (h2oliq > clwt)
then
155 onemrh= max( 1.e-10, 1.0-rhgrid )
156 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0)
158 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
159 tem2 = sqrt( sqrt(rhgrid) )
161 clouds1(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
173 qc_save(i,k) = qc(i,k)
174 qi_save(i,k) = qi(i,k)
175 qs_save(i,k) = qs(i,k)
179 if ( do_mynnedmf )
then
193 clouds1(i,k) = cldfra_bl(i,k)
196 if (qc(i,k) < 1.e-6 .and. cldfra_bl(i,k)>0.001)
then
200 if (nint(slmsk(i)) == 1)
then
201 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
203 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
207 clouds2(i,k) = max(0.0, qc(i,k) *
gfac * delp(i,k))
210 tc = t3d(i,k) - 273.15
213 snow_frac = min(0.5, max((p3d(i,k)-30000.0),0.0)/140000.0)
214 ice_frac = 1.0 - snow_frac
215 if (qi(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001)
then
216 qi(i,k) = ice_frac*qi_bl(i,k)
219 clouds5(i,k)=max(173.45 + 2.14*tc, 20.)
225 clouds4(i,k) = max(0.0, qi(i,k) *
gfac * delp(i,k))
228 if (qs(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001)
then
229 qs(i,k) = snow_frac*qi_bl(i,k)
232 clouds9(i,k)=max(2.*(173.45 + 2.14*tc), 50.)
235 clouds8(i,k) = max(0.0, qs(i,k) *
gfac * delp(i,k))
240 elseif (imp_physics /= imp_physics_gfdl)
then
249 if ( qi(i,k) > 1e-7 .OR. qc(i,k) > 1e-7 )
then
251 es = min( p3d(i,k), fpvs( t3d(i,k) ) )
252 qsat = max( qmin, eps * es / (p3d(i,k) + epsm1*es) )
253 rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
254 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k)
255 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
257 if (h2oliq > clwt)
then
258 onemrh= max( 1.e-10, 1.0-rhgrid )
259 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0)
261 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
262 tem2 = sqrt( sqrt(rhgrid) )
264 clouds1(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
276 if (imfdeepcnv == imfdeepcnv_gf .or. imfdeepcnv == imfdeepcnv_c3)
then
279 if ( qci_conv(i,k) > 0. )
then
284 liqfrac = min(1., max(0., (tk-244.)/29.))
285 qc(i,k) = qc(i,k)+qci_conv(i,k)*liqfrac
287 qi(i,k) = qi(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)
288 qs(i,k) = qs(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)
291 if (nint(slmsk(i)) == 1)
then
292 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
295 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
298 if(qi(i,k)>1.e-8)clouds5(i,k)=max( 173.45 + 2.14*tc , 20.)
299 if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*tc), 50.)
301 if ( conv_cf_opt .eq. 0 )
then
307 xl = xlv*liqfrac + xls*(1.-liqfrac)
308 tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
309 & - xlscp/exner(i,k)*qi(i,k)
311 es = min( p3d(i,k), fpvs( tlk ) )
312 qsat= max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
313 rsl = xl*qsat / (r_v*tlk**2)
315 qt = qc(i,k) + qi(i,k) + qv(i,k)
317 a = 1./(1. + xl*rsl/cpm)
324 sigq = 1.5e-3 * ud_mf(i,k)/dt * f
326 sigq = sqrt(sigq**2 + 1e-10)
327 qmq = a * (qt - qsat)
329 cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
330 if (qci_conv(i,k) .lt. 1e-9) cb_cf = 0.0
331 if (do_mynnedmf .and. qmq .ge. 0.0)
then
333 if (cb_cf .gt. 0.0)
then
334 clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
344 es = min( p3d(i,k), fpvs( t3d(i,k) ) )
345 qsat = max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
346 rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
347 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k)
348 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
350 if (h2oliq > clwt)
then
351 onemrh= max( 1.e-10, 1.0-rhgrid )
352 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0)
354 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
355 tem2 = sqrt( sqrt(rhgrid) )
357 clouds1(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
368 elseif (imfdeepcnv == imfdeepcnv_sas)
then
372 h2oliq = qlc(i,k)+qli(i,k)
373 if ( h2oliq > 0. )
then
378 liqfrac = min(1., max(0., (tk-244.)/29.))
380 qc(i,k) = qc(i,k)+qlc(i,k)
382 qi(i,k) = qi(i,k)+0.5*qli(i,k)
383 qs(i,k) = qs(i,k)+0.5*qli(i,k)
386 if (nint(slmsk(i)) == 1)
then
387 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
390 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
393 if(qi(i,k)>1.e-8)clouds5(i,k)=max( 173.45 + 2.14*tc , 20.)
394 if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*tc), 50.)
396 if ( conv_cf_opt .eq. 0 )
then
400 xl = xlv*liqfrac + xls*(1.-liqfrac)
401 tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
402 & - xlscp/exner(i,k)*qi(i,k)
404 es = min( p3d(i,k), fpvs( tlk ) )
405 qsat= max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
406 rsl = xl*qsat / (r_v*tlk**2)
408 qt = qc(i,k) + qi(i,k) + qv(i,k)
410 a = 1./(1. + xl*rsl/cpm)
417 sigq = 1.5e-3 * ud_mf(i,k)/dt * f
419 sigq = sqrt(sigq**2 + 1e-10)
420 qmq = a * (qt - qsat)
422 cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
423 if (h2oliq .lt. 1e-9) cb_cf = 0.0
424 if (do_mynnedmf .and. qmq .ge. 0.0)
then
426 if (cb_cf .gt. 0.0)
then
427 clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
437 es = min( p3d(i,k), fpvs( t3d(i,k) ) )
438 qsat = max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
439 rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
440 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k)
441 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
443 if (h2oliq > clwt)
then
444 onemrh= max( 1.e-10, 1.0-rhgrid )
445 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0)
447 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
448 tem2 = sqrt( sqrt(rhgrid) )
450 clouds1(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )