diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index 512257258..a7afa2ee0 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -81,7 +81,7 @@ end subroutine GFS_MP_generic_post_init !> \section gfs_mp_gen GFS MP Generic Post General Algorithm !> @{ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, & - imp_physics_thompson, imp_physics_mg, cal_pre, lssav, ldiag3d, cplflx, cplchm, con_g, dtf, frain, rainc, rain1, & + imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires, cal_pre, lssav, ldiag3d, cplflx, cplchm, con_g, dtf, frain, rainc, rain1, & rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_qv, rain0, ice0, snow0, & graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice, & totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, dt3dt, dq3dt, rain_cpl, rainc_cpl, snow_cpl, pwat, & @@ -93,7 +93,7 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt implicit none integer, intent(in) :: im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac - integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg + integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires logical, intent(in) :: cal_pre, lssav, ldiag3d, cplflx, cplchm real(kind=kind_phys), intent(in) :: dtf, frain, con_g @@ -179,6 +179,10 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt graupel = frain*graupel0 ! time-step graupel ice = frain*ice0 ! time-step ice snow = frain*snow0 ! time-step snow + + else if (imp_physics == imp_physics_fer_hires) then + tprcp = max (0.,rain) ! time-step convective and explicit precip + ice = frain*rain1*sr ! time-step ice end if if (lsm==lsm_ruc) then @@ -296,7 +300,6 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt ! if (snow0(i)+ice0(i)+graupel0(i)+csnow > 0.0) then ! Sfcprop%srflag(i) = 1. ! clu: set srflag to 'snow' (i.e. 1) ! endif -! compute fractional srflag total_precip = snow0(i)+ice0(i)+graupel0(i)+rain0(i)+rainc(i) if (total_precip > rainmin) then srflag(i) = (snow0(i)+ice0(i)+graupel0(i)+csnow)/total_precip diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 2e55b6ad5..3a11a9983 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -234,6 +234,14 @@ type = integer intent = in optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [cal_pre] standard_name = flag_for_precipitation_type_algorithm long_name = flag controls precip type algorithm diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index ec6134ed5..4bebae589 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -81,9 +81,9 @@ end subroutine GFS_PBL_generic_pre_finalize !! subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, & ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, & - ntwa, ntia, ntgl, ntoz, ntke, ntkev, trans_aero, ntchs, ntchm, & + ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm, & imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & - imp_physics_zhao_carr, imp_physics_mg, cplchm, ltaerosol, hybedmf, do_shoc, & + imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, cplchm, ltaerosol, hybedmf, do_shoc, & satmedmf, qgrs, vdftra, errmsg, errflg) use machine, only : kind_phys @@ -93,10 +93,10 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, integer, intent(in) :: im, levs, nvdiff, ntrac integer, intent(in) :: ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc - integer, intent(in) :: ntwa, ntia, ntgl, ntoz, ntke, ntkev, ntchs, ntchm + integer, intent(in) :: ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef,ntchs, ntchm logical, intent(in) :: trans_aero integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 - integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg + integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires logical, intent(in) :: cplchm, ltaerosol, hybedmf, do_shoc, satmedmf real(kind=kind_phys), dimension(im, levs, ntrac), intent(in) :: qgrs @@ -126,6 +126,20 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, vdftra(i,k,4) = qgrs(i,k,ntoz) enddo enddo + + ! Ferrier-Aligo + elseif (imp_physics == imp_physics_fer_hires) then + do k=1,levs + do i=1,im + vdftra(i,k,1) = qgrs(i,k,ntqv) + vdftra(i,k,2) = qgrs(i,k,ntcw) + vdftra(i,k,3) = qgrs(i,k,ntiw) + vdftra(i,k,4) = qgrs(i,k,ntrw) + vdftra(i,k,5) = qgrs(i,k,nqrimef) + vdftra(i,k,6) = qgrs(i,k,ntoz) + enddo + enddo + elseif (imp_physics == imp_physics_thompson) then ! Thompson if(ltaerosol) then @@ -263,9 +277,10 @@ end subroutine GFS_PBL_generic_post_finalize !! \htmlinclude GFS_PBL_generic_post_run.html !! subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, & - ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev, & + ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev,nqrimef, & trans_aero, ntchs, ntchm, & imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, imp_physics_mg, & + imp_physics_fer_hires, & ltaerosol, cplflx, cplchm, lssav, ldiag3d, lsidea, hybedmf, do_shoc, satmedmf, shinhong, do_ysu, & dvdftra, dusfc1, dvsfc1, dtsfc1, dqsfc1, dtf, dudt, dvdt, dtdt, htrsw, htrlw, xmu, & dqdt, dusfc_cpl, dvsfc_cpl, dtsfc_cpl, & @@ -280,10 +295,10 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, implicit none integer, intent(in) :: im, levs, nvdiff, ntrac, ntchs, ntchm - integer, intent(in) :: ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev + integer, intent(in) :: ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef logical, intent(in) :: trans_aero integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 - integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg + integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires logical, intent(in) :: ltaerosol, cplflx, cplchm, lssav, ldiag3d, lsidea logical, intent(in) :: hybedmf, do_shoc, satmedmf, shinhong, do_ysu @@ -365,6 +380,20 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqdt(i,k,ntoz) = dvdftra(i,k,4) enddo enddo + + elseif (imp_physics == imp_physics_fer_hires) then + ! Ferrier-Aligo + do k=1,levs + do i=1,im + dqdt(i,k,ntqv) = dvdftra(i,k,1) + dqdt(i,k,ntcw) = dvdftra(i,k,2) + dqdt(i,k,ntiw) = dvdftra(i,k,3) + dqdt(i,k,ntrw) = dvdftra(i,k,4) + dqdt(i,k,nqrimef) = dvdftra(i,k,5) + dqdt(i,k,ntoz) = dvdftra(i,k,6) + enddo + enddo + elseif (imp_physics == imp_physics_thompson) then ! Thompson if(ltaerosol) then diff --git a/physics/GFS_PBL_generic.meta b/physics/GFS_PBL_generic.meta index 25e696add..51764e04d 100644 --- a/physics/GFS_PBL_generic.meta +++ b/physics/GFS_PBL_generic.meta @@ -161,6 +161,14 @@ type = integer intent = in optional = F +[nqrimef] + standard_name = index_for_mass_weighted_rime_factor + long_name = tracer index for mass weighted rime factor + units = index + dimensions = () + type = integer + intent = in + optional = F [trans_aero] standard_name = flag_for_aerosol_convective_transport_and_PBL_diffusion long_name = flag for aerosol convective transport and PBL diffusion @@ -233,6 +241,14 @@ type = integer intent = in optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [cplchm] standard_name = flag_for_chemistry_coupling long_name = flag controlling cplchm collection (default off) @@ -473,6 +489,14 @@ type = integer intent = in optional = F +[nqrimef] + standard_name = index_for_mass_weighted_rime_factor + long_name = tracer index for mass weighted rime factor + units = index + dimensions = () + type = integer + intent = in + optional = F [trans_aero] standard_name = flag_for_aerosol_convective_transport_and_PBL_diffusion long_name = flag for aerosol convective transport and PBL diffusion @@ -545,6 +569,14 @@ type = integer intent = in optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [ltaerosol] standard_name = flag_for_aerosol_physics long_name = flag for aerosol physics diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 0f1dc291a..6098abefb 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -20,6 +20,7 @@ end subroutine GFS_rrtmg_pre_init subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd, Cldprop, Coupling, & Radtend, & ! input/output + f_ice, f_rain, f_rimef, flgmin, cwm, & ! F-A mp scheme only lm, im, lmk, lmp, & ! input kd, kt, kb, raddt, delp, dz, plvl, plyr, & ! output tlvl, tlyr, tsfg, tsfa, qlyr, olyr, & @@ -59,6 +60,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input & NSPC1 use module_radiation_clouds, only: NF_CLDS, & ! cld_init & progcld1, progcld3, & + & progcld2, & & progcld4, progcld5, & & progclduni use module_radsw_parameters, only: topfsw_type, sfcfsw_type, & @@ -80,8 +82,16 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input integer, intent(in) :: im, lm, lmk, lmp integer, intent(out) :: kd, kt, kb + +! F-A mp scheme only + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: f_ice + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: f_rain + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: f_rimef + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: cwm + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(in) :: flgmin real(kind=kind_phys), intent(out) :: raddt + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: delp real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: dz real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+1+LTP), intent(out) :: plvl @@ -520,7 +530,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ccnd(i,k,1) = tracer1(i,k,ntcw) ! liquid water/ice enddo enddo - elseif (Model%ncnd == 2) then ! MG + elseif (Model%ncnd == 2) then ! MG or F-A do k=1,LMK do i=1,IM ccnd(i,k,1) = tracer1(i,k,ntcw) ! liquid water @@ -749,6 +759,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Model%sup, Model%kdt, me, & clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs + elseif (Model%imp_physics == 11) then ! GFDL cloud scheme if (.not.Model%lgfdlmprad) then @@ -773,8 +784,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs endif - elseif(Model%imp_physics == 8 .or. Model%imp_physics == 6) then ! Thompson / WSM6 cloud micrphysics scheme - + elseif(Model%imp_physics == 8 .or. Model%imp_physics == 6 .or. & + Model%imp_physics == 15) then if (Model%kdt == 1) then Tbd%phy_f3d(:,:,Model%nleffr) = 10. Tbd%phy_f3d(:,:,Model%nieffr) = 50. @@ -795,7 +806,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! endif ! end_if_ntcw -! CCPP do k = 1, LMK do i = 1, IM clouds1(i,k) = clouds(i,k,1) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index d0c370882..7b40e2c1d 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -70,6 +70,51 @@ type = GFS_radtend_type intent = inout optional = F +[f_ice] + standard_name = fraction_of_ice_water_cloud + long_name = fraction of ice water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[f_rain] + standard_name = fraction_of_rain_water_cloud + long_name = fraction of rain water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[f_rimef] + standard_name = rime_factor + long_name = rime factor + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[flgmin] + standard_name = minimum_large_ice_fraction + long_name = minimum large ice fraction in F-A mp scheme + units = frac + dimensions = (2) + type = real + kind = kind_phys + intent = in + optional = F +[cwm] + standard_name = total_cloud_condensate_mixing_ratio_updated_by_physics + long_name = total cloud condensate mixing ratio (except water vapor) updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F [lm] standard_name = number_of_vertical_layers_for_radiation_calculations long_name = number of vertical layers for radiation calculation diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 6ec16f8b9..1e8545e98 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -14,21 +14,22 @@ end subroutine GFS_suite_interstitial_rad_reset_finalize !> \section arg_table_GFS_suite_interstitial_rad_reset_run Argument Table !! \htmlinclude GFS_suite_interstitial_rad_reset_run.html !! - subroutine GFS_suite_interstitial_rad_reset_run (Interstitial, errmsg, errflg) + subroutine GFS_suite_interstitial_rad_reset_run (Interstitial, Model, errmsg, errflg) - use GFS_typedefs, only: GFS_interstitial_type + use GFS_typedefs, only: GFS_control_type,GFS_interstitial_type implicit none ! interface variables type(GFS_interstitial_type), intent(inout) :: Interstitial + type(GFS_control_type), intent(in) :: Model character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg errmsg = '' errflg = 0 - call Interstitial%rad_reset() + call Interstitial%rad_reset(Model) end subroutine GFS_suite_interstitial_rad_reset_run @@ -459,11 +460,16 @@ end subroutine GFS_suite_interstitial_3_finalize !! \htmlinclude GFS_suite_interstitial_3_run.html !! #endif - subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & - ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, xlat, gq0, imp_physics, imp_physics_mg, imp_physics_zhao_carr,& - imp_physics_zhao_carr_pdf, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, prsi, prsl, prslk, rhcbot, & - rhcpbl, rhctop, rhcmax, islmsk, work1, work2, kpbl, kinver, & - clw, rhc, save_qc, save_qi, errmsg, errflg) + subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & + satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & + ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & + xlat, gq0, imp_physics, imp_physics_mg, & + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & + imp_physics_gfdl, imp_physics_thompson, & + imp_physics_wsm6, imp_physics_fer_hires, prsi, & + prsl, prslk, rhcbot,rhcpbl, rhctop, rhcmax, islmsk, & + work1, work2, kpbl, kinver,clw, rhc, save_qc, save_qi, & + errmsg, errflg) use machine, only: kind_phys @@ -472,7 +478,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, satmedmf, trans_tr ! interface variables integer, intent(in) :: im, levs, nn, ntrac, ntcw, ntiw, ntclamt, ntrw, & ntsw, ntrnc, ntsnc, ntgl, ntgnc, imp_physics, imp_physics_mg, imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & - imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 + imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires integer, dimension(im), intent(in) :: islmsk, kpbl, kinver logical, intent(in) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol @@ -619,7 +625,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, satmedmf, trans_tr else save_qi(:,:) = clw(:,:,1) endif - elseif (imp_physics == imp_physics_wsm6 .or. imp_physics == imp_physics_mg) then + elseif (imp_physics == imp_physics_wsm6 .or. imp_physics == imp_physics_mg .or. imp_physics == imp_physics_fer_hires) then do k=1,levs do i=1,im clw(i,k,1) = gq0(i,k,ntiw) ! ice @@ -656,7 +662,7 @@ end subroutine GFS_suite_interstitial_4_finalize subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, dtf, save_qc, save_qi, con_pi, & - gq0, clw, dqdti, errmsg, errflg) + gq0, clw, dqdti, imfdeepcnv, imfdeepcnv_gf, errmsg, errflg) use machine, only: kind_phys @@ -666,7 +672,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to integer, intent(in) :: im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, & ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imfdeepcnv, imfdeepcnv_gf logical, intent(in) :: ltaerosol, cplchm @@ -680,6 +686,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! dqdti may not be allocated real(kind=kind_phys), dimension(:,:), intent(inout) :: dqdti + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -729,7 +736,8 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to gq0(i,k,ntcw) = clw(i,k,2) ! water enddo enddo - if (imp_physics == imp_physics_thompson) then +! if (imp_physics == imp_physics_thompson) then + if (imp_physics == imp_physics_thompson .and. imfdeepcnv /= imfdeepcnv_gf) then if (ltaerosol) then do k=1,levs do i=1,im @@ -748,6 +756,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to enddo endif endif + else do k=1,levs do i=1,im diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index c07d9341a..e6e349a2a 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -9,6 +9,14 @@ type = GFS_interstitial_type intent = inout optional = F +[Model] + standard_name = GFS_control_type_instance + long_name = Fortran DDT containing FV3-GFS model control parameters + units = DDT + dimensions = () + type = GFS_control_type + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -1275,6 +1283,14 @@ type = integer intent = in optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [prsi] standard_name = air_pressure_at_interface long_name = air pressure at model layer interfaces @@ -1685,6 +1701,22 @@ kind = kind_phys intent = inout optional = F +[imfdeepcnv] + standard_name = flag_for_mass_flux_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imfdeepcnv_gf] + standard_name = flag_for_gf_deep_convection_scheme + long_name = flag for Grell-Freitas deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/cu_gf_deep.F90 b/physics/cu_gf_deep.F90 index 3e865c9ba..4afad80d1 100644 --- a/physics/cu_gf_deep.F90 +++ b/physics/cu_gf_deep.F90 @@ -14,7 +14,7 @@ module cu_gf_deep !> tuning constant for cloudwater/ice detrainment real(kind=kind_phys), parameter:: c1= 0.003 !.002 ! .0005 !> parameter to turn on or off evaporation of rainwater as done in sas - integer, parameter :: irainevap=0 + integer, parameter :: irainevap=1 !> max allowed fractional coverage (frh_thresh) real(kind=kind_phys), parameter::frh_thresh = .9 !> rh threshold. if fractional coverage ~ frh_thres, do not use cupa any further @@ -362,7 +362,7 @@ subroutine cu_gf_deep_run( & c1_max=c1 elocp=xlv/cp el2orc=xlv*xlv/(r_v*cp) - evfact=.2 + evfact=.4 ! .2 evfactl=.2 !evfact=.0 ! for 4F5f !evfactl=.4 @@ -1923,6 +1923,13 @@ subroutine cu_gf_deep_run( & ichoice,imid,ipr,itf,ktf, & its,ite, kts,kte, & dicycle,xf_dicycle ) + +!---------------evap below cloud base + + call rain_evap_below_cloudbase(itf,ktf,its,ite, & + kts,kte,ierr,kbcon,xmb,psur,xland,qo_cup, & + po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy) + k=1 do i=its,itf if(ierr(i).eq.0 .and.pre(i).gt.0.) then @@ -1971,7 +1978,7 @@ subroutine cu_gf_deep_run( & do k = ktop(i), 1, -1 rain = pwo(i,k) + edto(i) * pwdo(i,k) rn(i) = rn(i) + rain * xmb(i) * .001 * dtime - !if(po(i,k).gt.700.)then + if(po(i,k).gt.400.)then if(flg(i))then q1=qo(i,k)+(outq(i,k))*dtime t1=tn(i,k)+(outt(i,k))*dtime @@ -1996,7 +2003,7 @@ subroutine cu_gf_deep_run( & pre(i)=max(pre(i),0.) delqev(i) = delqev(i) + .001*dp*qevap(i)/g endif - !endif ! 700mb + endif ! 400mb endif enddo ! pre(i)=1000.*rn(i)/dtime @@ -2035,6 +2042,271 @@ end subroutine cu_gf_deep_run !> @} !>\ingroup cu_gf_deep_group + + + subroutine fct1d3 (ktop,n,dt,z,tracr,massflx,trflx_in,dellac,g) + +! --- modify a 1-D array of tracer fluxes for the purpose of maintaining +! --- monotonicity (including positive-definiteness) in the tracer field +! --- during tracer transport. + +! --- the underlying transport equation is (d tracr/dt) = - (d trflx/dz) +! --- where dz = |z(k+1)-z(k)| (k=1,...,n) and trflx = massflx * tracr +! --- physical dimensions of tracr,trflx,dz are arbitrary to some extent +! --- but are subject to the constraint dim[trflx] = dim[tracr*(dz/dt)]. + +! --- note: tracr is carried in grid cells while z and fluxes are carried on +! --- interfaces. interface variables at index k are at grid location k-1/2. +! --- sign convention: mass fluxes are considered positive in +k direction. + +! --- massflx and trflx_in must be provided independently to allow the +! --- algorithm to generate an auxiliary low-order (diffusive) tracer flux +! --- as a stepping stone toward the final product trflx_out. + + implicit none + integer,intent(in) :: n,ktop ! number of grid cells + real(kind=kind_phys) ,intent(in) :: dt,g ! transport time step + real(kind=kind_phys) ,intent(in) :: z(n+0) ! location of cell interfaces + real(kind=kind_phys) ,intent(in) :: tracr(n) ! the transported variable + real(kind=kind_phys) ,intent(in) :: massflx(n+0) ! mass flux across interfaces + real(kind=kind_phys) ,intent(in) :: trflx_in(n+0) ! original tracer flux + real(kind=kind_phys) ,intent(out):: dellac(n+0) ! modified tracr flux + real(kind=kind_phys) :: trflx_out(n+0) ! modified tracr flux + integer k,km1,kp1 + logical :: NaN, error=.false., vrbos=.true. + real(kind=kind_phys) dtovdz(n),trmax(n),trmin(n),flx_lo(n+0),antifx(n+0),clipped(n+0), & + soln_hi(n),totlin(n),totlout(n),soln_lo(n),clipin(n),clipout(n),arg + real(kind=kind_phys),parameter :: epsil=1.e-22 ! prevent division by zero + real(kind=kind_phys),parameter :: damp=1. ! damper of antidff flux (1=no damping) + NaN(arg) = .not. (arg.ge.0. .or. arg.le.0.) ! NaN detector + dtovdz(:)=0. + soln_lo(:)=0. + antifx(:)=0. + clipin(:)=0. + totlin(:)=0. + totlout(:)=0. + clipout(:)=0. + flx_lo(:)=0. + trmin(:)=0. + trmax(:)=0. + clipped(:)=0. + trflx_out(:)=0. + do k=1,ktop + dtovdz(k)=.01*dt/abs(z(k+1)-z(k))*g ! time step / grid spacing + if (z(k).eq.z(k+1)) error=.true. + end do +! if (vrbos .or. error) print '(a/(8es10.3))','(fct1d) dtovdz =',dtovdz + + do k=2,ktop + if (massflx(k).ge.0.) then + flx_lo(k)=massflx(k)*tracr(k-1) ! low-order flux, upstream + else + flx_lo(k)=massflx(k)*tracr(k) ! low-order flux, upstream + end if + antifx(k)=trflx_in(k)-flx_lo(k) ! antidiffusive flux + end do + flx_lo( 1)=trflx_in( 1) + flx_lo(ktop+1)=trflx_in(ktop+1) + antifx( 1)=0. + antifx(ktop+1)=0. +! --- clip low-ord fluxes to make sure they don't violate positive-definiteness + do k=1,ktop + totlout(k)=max(0.,flx_lo(k+1))-min(0.,flx_lo(k )) ! total flux out + clipout(k)=min(1.,tracr(k)/max(epsil,totlout(k))/ (1.0001*dtovdz(k))) + end do + + do k=2,ktop + if (massflx(k).ge.0.) then + flx_lo(k)=flx_lo(k)*clipout(k-1) + else + flx_lo(k)=flx_lo(k)*clipout(k) + end if + end do + if (massflx( 1).lt.0.) flx_lo( 1)=flx_lo( 1)*clipout(1) + if (massflx(ktop+1).gt.0.)flx_lo(ktop+1)=flx_lo(ktop+1)*clipout(ktop) + +! --- a positive-definite low-order (diffusive) solution can now be constructed + + do k=1,ktop + soln_lo(k)=tracr(k)-(flx_lo(k+1)-flx_lo(k))*dtovdz(k) ! low-ord solutn + dellac(k)=-(flx_lo(k+1)-flx_lo(k))*dtovdz(k)/dt + !dellac(k)=soln_lo(k) + end do + return + do k=1,ktop + km1=max(1,k-1) + kp1=min(ktop,k+1) + trmax(k)= max(soln_lo(km1),soln_lo(k),soln_lo(kp1), & + tracr (km1),tracr (k),tracr (kp1)) ! upper bound + trmin(k)=max(0.,min(soln_lo(km1),soln_lo(k),soln_lo(kp1), & + tracr (km1),tracr (k),tracr (kp1))) ! lower bound + end do + + do k=1,ktop + totlin (k)=max(0.,antifx(k ))-min(0.,antifx(k+1)) ! total flux in + totlout(k)=max(0.,antifx(k+1))-min(0.,antifx(k )) ! total flux out + + clipin (k)=min(damp,(trmax(k)-soln_lo(k))/max(epsil,totlin (k)) & + / (1.0001*dtovdz(k))) + clipout(k)=min(damp,(soln_lo(k)-trmin(k))/max(epsil,totlout(k)) & + / (1.0001*dtovdz(k))) + + if (NaN(clipin(k))) print *,'(fct1d) error: clipin is NaN, k=',k + if (NaN(clipout(k))) print *,'(fct1d) error: clipout is NaN, k=',k + + if (clipin(k).lt.0.) then +! print 100,'(fct1d) error: clipin < 0 at k =',k, & +! 'clipin',clipin(k),'trmax',trmax(k),'soln_lo',soln_lo(k), & +! 'totlin',totlin(k),'dt/dz',dtovdz(k) + error=.true. + end if + if (clipout(k).lt.0.) then +! print 100,'(fct1d) error: clipout < 0 at k =',k, & +! 'clipout',clipout(k),'trmin',trmin(k),'soln_lo',soln_lo(k), & +! 'totlout',totlout(k),'dt/dz',dtovdz(k) + error=.true. + end if +! 100 format (a,i3/(4(a10,"=",es9.2))) + end do + + do k=2,ktop + if (antifx(k).gt.0.) then + clipped(k)=antifx(k)*min(clipout(k-1),clipin(k)) + else + clipped(k)=antifx(k)*min(clipout(k),clipin(k-1)) + end if + trflx_out(k)=flx_lo(k)+clipped(k) + if (NaN(trflx_out(k))) then + print *,'(fct1d) error: trflx_out is NaN, k=',k + error=.true. + end if + end do + trflx_out( 1)=trflx_in( 1) + trflx_out(ktop+1)=trflx_in(ktop+1) + do k=1,ktop + soln_hi(k)=tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k) + dellac(k)=-g*(trflx_out(k+1)-trflx_out(k))*dtovdz(k)/dt + !dellac(k)=soln_hi(k) + end do + + if (vrbos .or. error) then +! do k=2,ktop +! write(32,99)k, & +! 'tracr(k)', tracr(k), & +! 'flx_in(k)', trflx_in(k), & +! 'flx_in(k+1)', trflx_in(k+1), & +! 'flx_lo(k)', flx_lo(k), & +! 'flx_lo(k+1)', flx_lo(k+1), & +! 'soln_lo(k)', soln_lo(k), & +! 'trmin(k)', trmin(k), & +! 'trmax(k)', trmax(k), & +! 'totlin(k)', totlin(k), & +! 'totlout(k)', totlout(k), & +! 'clipin(k-1)', clipin(k-1), & +! 'clipin(k)', clipin(k), & +! 'clipout(k-1)', clipout(k-1), & +! 'clipout(k)', clipout(k), & +! 'antifx(k)', antifx(k), & +! 'antifx(k+1)', antifx(k+1), & +! 'clipped(k)', clipped(k), & +! 'clipped(k+1)', clipped(k+1), & +! 'flx_out(k)', trflx_out(k), & +! 'flx_out(k+1)', trflx_out(k+1), & +! 'dt/dz(k)', dtovdz(k), & +! 'final', tracr(k)-(trflx_out(k+1)-trflx_out(k))*dtovdz(k) +! 99 format ('(trc1d) k =',i4/(3(a13,'=',es13.6))) +! end do + if (error) stop '(fct1d error)' + end if + + return + end subroutine fct1d3 + + subroutine rain_evap_below_cloudbase(itf,ktf, its,ite, kts,kte,ierr, & + kbcon,xmb,psur,xland,qo_cup, & + po_cup,qes_cup,pwavo,edto,pwevo,pre,outt,outq) !,outbuoy) + + implicit none + real(kind=kind_phys), parameter :: alp1=5.44e-4 & !1/sec + ,alp2=5.09e-3 & !unitless + ,alp3=0.5777 & !unitless + ,c_conv=0.05 !conv fraction area, unitless + + + integer ,intent(in) :: itf,ktf, its,ite, kts,kte + integer, dimension(its:ite) ,intent(in) :: ierr,kbcon + real(kind=kind_phys), dimension(its:ite) ,intent(in) ::psur,xland,pwavo,edto,pwevo,xmb + real(kind=kind_phys), dimension(its:ite,kts:kte),intent(in) :: po_cup,qo_cup,qes_cup + real(kind=kind_phys), dimension(its:ite) ,intent(inout) :: pre + real(kind=kind_phys), dimension(its:ite,kts:kte),intent(inout) :: outt,outq !,outbuoy + + !real, dimension(its:ite) ,intent(out) :: tot_evap_bcb + !real, dimension(its:ite,kts:kte),intent(out) :: evap_bcb,net_prec_bcb + + !-- locals + integer :: i,k + real(kind=kind_phys) :: RH_cr , del_t,del_q,dp,q_deficit + real(kind=kind_phys), dimension(its:ite,kts:kte) :: evap_bcb,net_prec_bcb + real(kind=kind_phys), dimension(its:ite) :: tot_evap_bcb + + do i=its,itf + evap_bcb (i,:)= 0.0 + net_prec_bcb(i,:)= 0.0 + tot_evap_bcb(i) = 0.0 + if(ierr(i) /= 0) cycle + + !-- critical rel humidity + RH_cr=0.9*xland(i)+0.7*(1-xland(i)) + !RH_cr=1. + + !-- net precipitation (after downdraft evap) at cloud base, available to + !evap + k=kbcon(i) + !net_prec_bcb(i,k) = xmb(i)*(pwavo(i)+edto(i)*pwevo(i)) !-- pwevo<0. + net_prec_bcb(i,k) = pre(i) + + do k=kbcon(i)-1, kts, -1 + + q_deficit = max(0.,(RH_cr*qes_cup(i,k) -qo_cup(i,k))) + + if(q_deficit < 1.e-6) then + net_prec_bcb(i,k)= net_prec_bcb(i,k+1) + cycle + endif + + dp = 100.*(po_cup(i,k)-po_cup(i,k+1)) + + !--units here: kg[water]/kg[air}/sec + evap_bcb(i,k) = c_conv * alp1 * q_deficit * & + ( sqrt(po_cup(i,k)/psur(i))/alp2 *net_prec_bcb(i,k+1)/c_conv )**alp3 + + !--units here: kg[water]/kg[air}/sec * kg[air]/m3 * m = kg[water]/m2/sec + evap_bcb(i,k)= evap_bcb(i,k)*dp/g + + if((net_prec_bcb(i,k+1) - evap_bcb(i,k)).lt.0.) cycle + if((pre(i) - evap_bcb(i,k)).lt.0.) cycle + net_prec_bcb(i,k)= net_prec_bcb(i,k+1) - evap_bcb(i,k) + + tot_evap_bcb(i) = tot_evap_bcb(i)+evap_bcb(i,k) + + !-- feedback + del_q = evap_bcb(i,k)*g/dp ! > 0., units: kg[water]/kg[air}/sec + del_t = -evap_bcb(i,k)*g/dp*(xlv/cp) ! < 0., units: K/sec + +! print*,"ebcb2",k,del_q*86400,del_t*86400 + + outq (i,k) = outq (i,k) + del_q + outt (i,k) = outt (i,k) + del_t + !outbuoy(i,k) = outbuoy(i,k) + cp*del_t+xlv*del_q + + pre(i) = pre(i) - evap_bcb(i,k) + enddo + enddo + + end subroutine rain_evap_below_cloudbase + + + subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & pw,ccn,pwev,edtmax,edtmin,edtc,psum2,psumh, & rho,aeroevap,itf,ktf, & @@ -2747,9 +3019,8 @@ subroutine cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2 xff_ens3(12)=0. xff_ens3(13)= 0. xff_ens3(16)= 0. -! closure_n(i)=12. -! hli 05/01/2018 closure_n(i)=12. -! xff_dicycle = 0. +! closure_n(i)=12. +! xff_dicycle = 0. endif !xff0 endif ! ichoice @@ -3682,7 +3953,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & prop_b(kts:kte)=0 iall=0 c0=.002 - clwdet=100. + clwdet=50. bdsp=bdispm ! !--- no precip for small clouds diff --git a/physics/cu_gf_driver.F90 b/physics/cu_gf_driver.F90 index 58a30749a..53e26fb46 100644 --- a/physics/cu_gf_driver.F90 +++ b/physics/cu_gf_driver.F90 @@ -7,8 +7,9 @@ module cu_gf_driver ! DH* TODO: replace constants with arguments to cu_gf_driver_run use physcons , g => con_g, cp => con_cp, xlv => con_hvap, r_v => con_rv use machine , only: kind_phys - use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap + use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3 use cu_gf_sh , only: cu_gf_sh_run + use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber implicit none @@ -68,11 +69,12 @@ end subroutine cu_gf_driver_finalize !! !>\section gen_gf_driver GSD GF Cumulus Scheme General Algorithm !> @{ - subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & - forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d, & - us,vs,t2di,w,qv2di_spechum,p2di,psuri, & - hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & - pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & + subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & + forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d, & + us,vs,t2di,w,qv2di_spechum,p2di,psuri, & + hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & + pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & + nwfa,con_rd,gq0,ntinc,ntlnc,imp_physics,imp_physics_thompson, & errmsg,errflg) !------------------------------------------------------------- implicit none @@ -94,7 +96,7 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & integer :: ishallow_g3 ! depend on imfshalcnv !------------------------------------------------------------- integer :: its,ite, jts,jte, kts,kte - integer, intent(in ) :: im,ix,km + integer, intent(in ) :: im,ix,km,ntracer real(kind=kind_phys), dimension( ix , km ), intent(in ) :: forcet,forceqv_spechum,w,phil real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: t,us,vs @@ -104,16 +106,16 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & real(kind=kind_phys), dimension( ix , km ), intent(out ) :: cnvw_moist,cnvc real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: cliw, clcw -!hj change from ix to im +! change from ix to im integer, dimension (im), intent(inout) :: hbot,htop,kcnv integer, dimension (im), intent(in) :: xland real(kind=kind_phys), dimension (im), intent(in) :: pbl integer, dimension (ix) :: tropics -! ruc variable +! ruc variable real(kind=kind_phys), dimension (im) :: hfx2,qfx2,psuri real(kind=kind_phys), dimension (im,km) :: ud_mf,dd_mf,dt_mf real(kind=kind_phys), dimension (im), intent(inout) :: raincv,cld1d -!hj end change ix to im +! end change ix to im real(kind=kind_phys), dimension (ix,km) :: t2di,p2di ! Specific humidity from FV3 real(kind=kind_phys), dimension (ix,km), intent(in) :: qv2di_spechum @@ -123,80 +125,76 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! real(kind=kind_phys), dimension( im ),intent(in) :: garea real(kind=kind_phys), intent(in ) :: dt + +! additional variables for number concentrations + real(kind=kind_phys), intent(in) :: nwfa(1:im,1:km) + real(kind=kind_phys), intent(in) :: con_rd + real(kind=kind_phys), dimension(im,km,ntracer), intent(inout) :: gq0 + integer, intent(in) :: imp_physics,imp_physics_thompson,ntlnc,ntinc + integer, intent(in ) :: imfshalcnv character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg -!hj define locally for now. - integer, dimension(im),intent(inout) :: cactiv ! hli for gf -!hj change from ix to im +! define locally for now. + integer, dimension(im),intent(inout) :: cactiv integer, dimension(im) :: k22_shallow,kbcon_shallow,ktop_shallow real(kind=kind_phys), dimension(im) :: ht -!hj change -! -!+lxz -!hj real(kind=kind_phys) :: dx real(kind=kind_phys), dimension(im) :: dx -! local vars -!hj change ix to im - real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws - real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm - real(kind=kind_phys), dimension (im,km) :: outts,outqs,outqcs,outu,outv,outus,outvs - real(kind=kind_phys), dimension (im,km) :: outtm,outqm,outqcm,submm,cupclwm - real(kind=kind_phys), dimension (im,km) :: cnvwt,cnvwts,cnvwtm - real(kind=kind_phys), dimension (im,km) :: hco,hcdo,zdo,zdd,hcom,hcdom,zdom - real(kind=kind_phys), dimension (km) :: zh - real(kind=kind_phys), dimension (im) :: tau_ecmwf,edt,edtm,edtd,ter11,aa0,xlandi - real(kind=kind_phys), dimension (im) :: pret,prets,pretm,hexec - real(kind=kind_phys), dimension (im,10) :: forcing,forcing2 -!+lxz - integer, dimension (im) :: kbcon, ktop,ierr,ierrs,ierrm,kpbli - integer, dimension (im) :: k22s,kbcons,ktops,k22,jmin,jminm - integer, dimension (im) :: kbconm,ktopm,k22m -!hj end change ix to im -!.lxz - integer :: iens,ibeg,iend,jbeg,jend,n - integer :: ibegh,iendh,jbegh,jendh - integer :: ibegc,iendc,jbegc,jendc,kstop - real(kind=kind_phys) :: rho_dryar,temp - real(kind=kind_phys) :: pten,pqen,paph,zrho,pahfs,pqhfl,zkhvfl,pgeoh -!hj 10/11/2016: ipn is an input in fim. set it to zero here. - integer, parameter :: ipn = 0 + real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws + real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm + real(kind=kind_phys), dimension (im,km) :: outts,outqs,outqcs,outu,outv,outus,outvs + real(kind=kind_phys), dimension (im,km) :: outtm,outqm,outqcm,submm,cupclwm + real(kind=kind_phys), dimension (im,km) :: cnvwt,cnvwts,cnvwtm + real(kind=kind_phys), dimension (im,km) :: hco,hcdo,zdo,zdd,hcom,hcdom,zdom + real(kind=kind_phys), dimension (km) :: zh + real(kind=kind_phys), dimension (im) :: tau_ecmwf,edt,edtm,edtd,ter11,aa0,xlandi + real(kind=kind_phys), dimension (im) :: pret,prets,pretm,hexec + real(kind=kind_phys), dimension (im,10) :: forcing,forcing2 + + integer, dimension (im) :: kbcon, ktop,ierr,ierrs,ierrm,kpbli + integer, dimension (im) :: k22s,kbcons,ktops,k22,jmin,jminm + integer, dimension (im) :: kbconm,ktopm,k22m + + integer :: iens,ibeg,iend,jbeg,jend,n + integer :: ibegh,iendh,jbegh,jendh + integer :: ibegc,iendc,jbegc,jendc,kstop + real(kind=kind_phys), dimension(im,km) :: rho_dryar + real(kind=kind_phys) :: pten,pqen,paph,zrho,pahfs,pqhfl,zkhvfl,pgeoh + integer, parameter :: ipn = 0 ! ! basic environmental input includes moisture convergence (mconv) ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off ! convection for this call only and at that particular gridpoint ! -!hj 10/11/2016: change ix to im. - real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi - real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg - real(kind=kind_phys), dimension (im) :: ccn,z1,psur,cuten,cutens,cutenm - real(kind=kind_phys), dimension (im) :: umean,vmean,pmean - real(kind=kind_phys), dimension (im) :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv -!hj end change ix to im - - integer :: i,j,k,icldck,ipr,jpr,jpr_deep,ipr_deep - integer :: itf,jtf,ktf,iss,jss,nbegin,nend - integer :: high_resolution - real(kind=kind_phys) :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter - real(kind=kind_phys) :: dsubclw,dsubclws,dsubclwm,ztm,ztq,hfm,qfm,rkbcon,rktop !-lxz -!hj change ix to im - real(kind=kind_phys), dimension (im) :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep - character*50 :: ierrc(im),ierrcm(im) - character*50 :: ierrcs(im) -!hj end change ix to im -! ruc variable -!hj hfx2 -- sensible heat flux (k m/s), positive upward from sfc -!hj qfx2 -- latent heat flux (kg/kg m/s), positive upward from sfc -!hj gf needs them in w/m2. define hfx and qfx after simple unit conversion - real(kind=kind_phys), dimension (im) :: hfx,qfx - real(kind=kind_phys) tem,tem1,tf,tcr,tcrf - - parameter (tf=243.16, tcr=270.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=263.16, tcr=273.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - !parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf)) ! as fim - ! initialize ccpp error handling variables + real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi + real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg + real(kind=kind_phys), dimension (im) :: ccn,z1,psur,cuten,cutens,cutenm + real(kind=kind_phys), dimension (im) :: umean,vmean,pmean + real(kind=kind_phys), dimension (im) :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv + + integer :: i,j,k,icldck,ipr,jpr,jpr_deep,ipr_deep + integer :: itf,jtf,ktf,iss,jss,nbegin,nend + integer :: high_resolution + real(kind=kind_phys) :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter + real(kind=kind_phys) :: dsubclw,dsubclws,dsubclwm,dtime_max,ztm,ztq,hfm,qfm,rkbcon,rktop + real(kind=kind_phys), dimension(km) :: massflx,trcflx_in1,clw_in1,clw_ten1,po_cup +! real(kind=kind_phys), dimension(km) :: trcflx_in2,clw_in2,clw_ten2 + real(kind=kind_phys), dimension (im) :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep + character*50 :: ierrc(im),ierrcm(im) + character*50 :: ierrcs(im) +! ruc variable +! hfx2 -- sensible heat flux (k m/s), positive upward from sfc +! qfx2 -- latent heat flux (kg/kg m/s), positive upward from sfc +! gf needs them in w/m2. define hfx and qfx after simple unit conversion + real(kind=kind_phys), dimension (im) :: hfx,qfx + real(kind=kind_phys) tem,tem1,tf,tcr,tcrf + + parameter (tf=243.16, tcr=270.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=263.16, tcr=273.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) + !parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf)) ! as fim + ! initialize ccpp error handling variables errmsg = '' errflg = 0 ! @@ -212,8 +210,7 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! ! these should be coming in from outside ! -! print*,'hli in gf cactiv',cactiv -! cactiv(:) = 0 +! cactiv(:) = 0 rand_mom(:) = 0. rand_vmas(:) = 0. rand_clos(:,:) = 0. @@ -232,112 +229,113 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! !> - Set tuning constants for radiation coupling ! - tun_rad_shall(:)=.02 - tun_rad_mid(:)=.15 - tun_rad_deep(:)=.13 - edt(:)=0. - edtm(:)=0. - edtd(:)=0. - zdd(:,:)=0. - flux_tun(:)=5. -!hj 10/11/2016 dx and tscl_kf are replaced with input dx(i), is dlength. - ! dx for scale awareness -!hj dx=40075000./float(lonf) -!hj tscl_kf=dx/25000. - ccn(its:ite)=150. - ! - if (imfshalcnv == 3) then - ishallow_g3 = 1 - else - ishallow_g3 = 0 - end if - high_resolution=0 - subcenter=0. - iens=1 + tun_rad_shall(:)=.02 + tun_rad_mid(:)=.15 + tun_rad_deep(:)=.13 + edt(:)=0. + edtm(:)=0. + edtd(:)=0. + zdd(:,:)=0. + flux_tun(:)=5. +! 10/11/2016 dx and tscl_kf are replaced with input dx(i), is dlength. +! dx for scale awareness +! dx=40075000./float(lonf) +! tscl_kf=dx/25000. + ccn(its:ite)=150. + + if (imfshalcnv == 3) then + ishallow_g3 = 1 + else + ishallow_g3 = 0 + end if + high_resolution=0 + subcenter=0. + iens=1 ! ! these can be set for debugging ! - ipr=0 - jpr=0 - ipr_deep=0 - jpr_deep= 0 !53322 ! 528196 !0 ! 1136 !0 !421755 !3536 + ipr=0 + jpr=0 + ipr_deep=0 + jpr_deep= 0 !53322 ! 528196 !0 ! 1136 !0 !421755 !3536 ! ! - ibeg=its - iend=ite - tcrit=258. - - ztm=0. - ztq=0. - hfm=0. - qfm=0. - ud_mf =0. - dd_mf =0. - dt_mf =0. - tau_ecmwf(:)=0. + ibeg=its + iend=ite + tcrit=258. + + ztm=0. + ztq=0. + hfm=0. + qfm=0. + ud_mf =0. + dd_mf =0. + dt_mf =0. + tau_ecmwf(:)=0. ! - j=1 - ht(:)=phil(:,1)/g - do i=its,ite - cld1d(i)=0. - zo(i,:)=phil(i,:)/g - dz8w(i,1)=zo(i,2)-zo(i,1) - zh(1)=0. - kpbli(i)=2 - do k=kts+1,ktf - dz8w(i,k)=zo(i,k+1)-zo(i,k) - enddo - do k=kts+1,ktf - zh(k)=zh(k-1)+dz8w(i,k-1) - if(zh(k).gt.pbl(i))then - kpbli(i)=max(2,k) - exit - endif - enddo - enddo - do i= its,itf - forcing(i,:)=0. - forcing2(i,:)=0. - ccn(i)=100. - hbot(i) =kte - htop(i) =kts - raincv(i)=0. - xlandi(i)=real(xland(i)) -! if(abs(xlandi(i)-1.).le.1.e-3) tun_rad_shall(i)=.15 -! if(abs(xlandi(i)-1.).le.1.e-3) flux_tun(i)=1.5 + j=1 + ht(:)=phil(:,1)/g + do i=its,ite + cld1d(i)=0. + zo(i,:)=phil(i,:)/g + dz8w(i,1)=zo(i,2)-zo(i,1) + zh(1)=0. + kpbli(i)=2 + do k=kts+1,ktf + dz8w(i,k)=zo(i,k+1)-zo(i,k) + enddo + do k=kts+1,ktf + zh(k)=zh(k-1)+dz8w(i,k-1) + if(zh(k).gt.pbl(i))then + kpbli(i)=max(2,k) + exit + endif + enddo enddo + do i= its,itf - mconv(i)=0. + forcing(i,:)=0. + forcing2(i,:)=0. + ccn(i)=100. + hbot(i) =kte + htop(i) =kts + raincv(i)=0. + xlandi(i)=real(xland(i)) +! if(abs(xlandi(i)-1.).le.1.e-3) tun_rad_shall(i)=.15 +! if(abs(xlandi(i)-1.).le.1.e-3) flux_tun(i)=1.5 enddo - do k=kts,kte do i= its,itf - omeg(i,k)=0. - zu(i,k)=0. - zum(i,k)=0. - zus(i,k)=0. - zd(i,k)=0. - zdm(i,k)=0. + mconv(i)=0. enddo + do k=kts,kte + do i= its,itf + omeg(i,k)=0. + zu(i,k)=0. + zum(i,k)=0. + zus(i,k)=0. + zd(i,k)=0. + zdm(i,k)=0. + enddo enddo psur(:)=0.01*psuri(:) do i=its,itf - ter11(i)=max(0.,ht(i)) + ter11(i)=max(0.,ht(i)) enddo do k=kts,kte - do i=its,ite - cnvw(i,k)=0. - cnvc(i,k)=0. - gdc(i,k,1)=0. - gdc(i,k,2)=0. - gdc(i,k,3)=0. - gdc(i,k,4)=0. - gdc(i,k,7)=0. - gdc(i,k,8)=0. - gdc(i,k,9)=0. - gdc(i,k,10)=0. - gdc2(i,k,1)=0. - enddo + do i=its,ite + cnvw(i,k)=0. + cnvc(i,k)=0. + gdc(i,k,1)=0. + gdc(i,k,2)=0. + gdc(i,k,3)=0. + gdc(i,k,4)=0. + gdc(i,k,7)=0. + gdc(i,k,8)=0. + gdc(i,k,9)=0. + gdc(i,k,10)=0. + gdc2(i,k,1)=0. + enddo enddo ierr(:)=0 ierrm(:)=0 @@ -410,88 +408,80 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & subm(:,:)=0. dhdt(:,:)=0. - !print*,'hli t2di',t2di - !print*,'hli forcet',forcet do k=kts,ktf - do i=its,itf - p2d(i,k)=0.01*p2di(i,k) - po(i,k)=p2d(i,k) !*.01 - rhoi(i,k) = 100.*p2d(i,k)/(287.04*(t2di(i,k)*(1.+0.608*qv2di(i,k)))) - qcheck(i,k)=qv(i,k) - tn(i,k)=t(i,k)!+forcet(i,k)*dt - qo(i,k)=max(1.e-16,qv(i,k))!+forceqv(i,k)*dt - t2d(i,k)=t2di(i,k)-forcet(i,k)*dt - !print*,'hli t2di(i,k),forcet(i,k),dt,t2d(i,k)',t2di(i,k),forcet(i,k),dt,t2d(i,k) - q2d(i,k)=max(1.e-16,qv2di(i,k)-forceqv(i,k)*dt) - if(qo(i,k).lt.1.e-16)qo(i,k)=1.e-16 - tshall(i,k)=t2d(i,k) - qshall(i,k)=q2d(i,k) -!hj if(ipn.eq.jpr_deep)then -!hj write(12,123)k,dt,p2d(i,k),t2d(i,k),tn(i,k),q2d(i,k),qo(i,k),forcet(i,k) -!hj endif - enddo + do i=its,itf + p2d(i,k)=0.01*p2di(i,k) + po(i,k)=p2d(i,k) !*.01 + rhoi(i,k) = 100.*p2d(i,k)/(287.04*(t2di(i,k)*(1.+0.608*qv2di(i,k)))) + qcheck(i,k)=qv(i,k) + tn(i,k)=t(i,k)!+forcet(i,k)*dt + qo(i,k)=max(1.e-16,qv(i,k))!+forceqv(i,k)*dt + t2d(i,k)=t2di(i,k)-forcet(i,k)*dt + q2d(i,k)=max(1.e-16,qv2di(i,k)-forceqv(i,k)*dt) + if(qo(i,k).lt.1.e-16)qo(i,k)=1.e-16 + tshall(i,k)=t2d(i,k) + qshall(i,k)=q2d(i,k) + enddo enddo 123 format(1x,i2,1x,2(1x,f8.0),1x,2(1x,f8.3),3(1x,e13.5)) do i=its,itf - do k=kts,kpbli(i) + do k=kts,kpbli(i) tshall(i,k)=t(i,k) qshall(i,k)=max(1.e-16,qv(i,k)) - enddo + enddo enddo ! -!hj converting hfx2 and qfx2 to w/m2 -!hj hfx=cp*rho*hfx2 -!hj qfx=xlv*qfx2 +! converting hfx2 and qfx2 to w/m2 +! hfx=cp*rho*hfx2 +! qfx=xlv*qfx2 do i=its,itf - hfx(i)=hfx2(i)*cp*rhoi(i,1) - qfx(i)=qfx2(i)*xlv*rhoi(i,1) - dx(i) = sqrt(garea(i)) - !print*,'hli dx', dx(i) + hfx(i)=hfx2(i)*cp*rhoi(i,1) + qfx(i)=qfx2(i)*xlv*rhoi(i,1) + dx(i) = sqrt(garea(i)) enddo -!hj write(0,*),'hfx',hfx(3),qfx(3),rhoi(3,1) -!hj + do i=its,itf - do k=kts,kpbli(i) - tn(i,k)=t(i,k) - qo(i,k)=max(1.e-16,qv(i,k)) - enddo + do k=kts,kpbli(i) + tn(i,k)=t(i,k) + qo(i,k)=max(1.e-16,qv(i,k)) + enddo enddo nbegin=0 nend=0 - do i=its,itf - do k=kts,kpbli(i) - dhdt(i,k)=cp*(forcet(i,k)+(t(i,k)-t2di(i,k))/dt) + & - xlv*(forceqv(i,k)+(qv(i,k)-qv2di(i,k))/dt) -! tshall(i,k)=t(i,k) -! qshall(i,k)=qv(i,k) - enddo - enddo - do k= kts+1,ktf-1 - do i = its,itf - if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then - dp=-.5*(p2d(i,k+1)-p2d(i,k-1)) - umean(i)=umean(i)+us(i,k)*dp - vmean(i)=vmean(i)+vs(i,k)*dp - pmean(i)=pmean(i)+dp - endif - enddo + do i=its,itf + do k=kts,kpbli(i) + dhdt(i,k)=cp*(forcet(i,k)+(t(i,k)-t2di(i,k))/dt) + & + xlv*(forceqv(i,k)+(qv(i,k)-qv2di(i,k))/dt) +! tshall(i,k)=t(i,k) +! qshall(i,k)=qv(i,k) enddo - do k=kts,ktf-1 + enddo + do k= kts+1,ktf-1 do i = its,itf - omeg(i,k)= w(i,k) !-g*rhoi(i,k)*w(i,k) -! dq=(q2d(i,k+1)-q2d(i,k)) -! mconv(i)=mconv(i)+omeg(i,k)*dq/g - enddo + if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then + dp=-.5*(p2d(i,k+1)-p2d(i,k-1)) + umean(i)=umean(i)+us(i,k)*dp + vmean(i)=vmean(i)+vs(i,k)*dp + pmean(i)=pmean(i)+dp + endif enddo + enddo + do k=kts,ktf-1 do i = its,itf - if(mconv(i).lt.0.)mconv(i)=0. + omeg(i,k)= w(i,k) !-g*rhoi(i,k)*w(i,k) +! dq=(q2d(i,k+1)-q2d(i,k)) +! mconv(i)=mconv(i)+omeg(i,k)*dq/g enddo + enddo + do i = its,itf + if(mconv(i).lt.0.)mconv(i)=0. + enddo ! !---- call cumulus parameterization ! if(ishallow_g3.eq.1)then -! + do i=its,ite ierrs(i)=0 ierrm(i)=0 @@ -499,14 +489,13 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ! !> - Call shallow: cu_gf_sh_run() ! - ! print*,'hli bf shallow t2d',t2d call cu_gf_sh_run (us,vs, & ! input variables, must be supplied zo,t2d,q2d,ter11,tshall,qshall,p2d,psur,dhdt,kpbli, & - rhoi,hfx,qfx,xlandi,ichoice_s,tcrit,dt, & + rhoi,hfx,qfx,xlandi,ichoice_s,tcrit,dt, & ! input variables. ierr should be initialized to zero or larger than zero for ! turning off shallow convection for grid points - zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs, & + zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs, & ! output tendencies outts,outqs,outqcs,outus,outvs,cnvwt,prets,cupclws, & ! dimesnional variables @@ -524,8 +513,8 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ipr=0 jpr_deep=0 !340765 !> - Call cu_gf_deep_run() for middle GF convection - if(imid_gf == 1)then - call cu_gf_deep_run( & + if(imid_gf == 1)then + call cu_gf_deep_run( & itf,ktf,its,ite, kts,kte & ,dicycle_m & ,ichoicem & @@ -594,16 +583,16 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & ,jminm,tropics) do i=its,itf - do k=kts,ktf + do k=kts,ktf qcheck(i,k)=qv(i,k) +outqs(i,k)*dt - enddo + enddo enddo !> - Call neg_check() for middle GF convection call neg_check('mid',ipn,dt,qcheck,outqm,outtm,outum,outvm, & outqcm,pretm,its,ite,kts,kte,itf,ktf,ktopm) - endif + endif !> - Call cu_gf_deep_run() for deep GF convection - if(ideep.eq.1)then + if(ideep.eq.1)then call cu_gf_deep_run( & itf,ktf,its,ite, kts,kte & @@ -673,15 +662,15 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & #endif ,k22 & ,jmin,tropics) - jpr=0 - ipr=0 - do i=its,itf - do k=kts,ktf - qcheck(i,k)=qv(i,k) +(outqs(i,k)+outqm(i,k))*dt - enddo - enddo + jpr=0 + ipr=0 + do i=its,itf + do k=kts,ktf + qcheck(i,k)=qv(i,k) +(outqs(i,k)+outqm(i,k))*dt + enddo + enddo !> - Call neg_check() for deep GF convection - call neg_check('deep',ipn,dt,qcheck,outq,outt,outu,outv, & + call neg_check('deep',ipn,dt,qcheck,outq,outt,outu,outv, & outqc,pret,its,ite,kts,kte,itf,ktf,ktop) ! endif @@ -730,6 +719,11 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & enddo ! do i=its,itf + massflx(:)=0. + trcflx_in1(:)=0. + clw_in1(:)=0. + clw_ten1(:)=0. + po_cup(:)=0. kstop=kts if(ktopm(i).gt.kts .or. ktop(i).gt.kts)kstop=max(ktopm(i),ktop(i)) if(ktops(i).gt.kts)kstop=max(kstop,ktops(i)) @@ -738,7 +732,8 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & if(kbcon(i).gt.2 .or. kbconm(i).gt.2)then hbot(i)=max(kbconm(i),kbcon(i)) !jmin(i) endif -!kbcon(i) + + dtime_max=dt do k=kts,kstop cnvc(i,k) = 0.04 * log(1. + 675. * zu(i,k) * xmb(i)) + & 0.04 * log(1. + 675. * zum(i,k) * xmbm(i)) + & @@ -754,66 +749,117 @@ subroutine cu_gf_driver_run(garea,im,ix,km,dt,cactiv, & us(i,k)=us(i,k)+outu(i,k)*cuten(i)*dt +outum(i,k)*cutenm(i)*dt +outus(i,k)*cutens(i)*dt vs(i,k)=vs(i,k)+outv(i,k)*cuten(i)*dt +outvm(i,k)*cutenm(i)*dt +outvs(i,k)*cutens(i)*dt -!hj 10/11/2016: don't need gdc and gdc2 yet for gsm. -!hli 08/18/2017: couple gdc to radiation - gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i)) ! my mod + gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i)) ! my mod gdc2(i,k,1)=max(0.,tun_rad_deep(i)*(cupclwm(i,k)*cutenm(i)+cupclw(i,k)*cuten(i))) gdc(i,k,2)=(outt(i,k))*86400. gdc(i,k,3)=(outtm(i,k))*86400. gdc(i,k,4)=(outts(i,k))*86400. gdc(i,k,7)=-(gdc(i,k,7)-sqrt(us(i,k)**2 +vs(i,k)**2))/dt - !gdc(i,k,8)=(outq(i,k))*86400.*xlv/cp + !gdc(i,k,8)=(outq(i,k))*86400.*xlv/cp gdc(i,k,8)=(outqm(i,k)+outqs(i,k)+outq(i,k))*86400.*xlv/cp gdc(i,k,9)=gdc(i,k,2)+gdc(i,k,3)+gdc(i,k,4) - if((gdc(i,k,1).ge.0.5).or.(gdc2(i,k,1).ge.0.5))then - print*,'hli gdc(i,k,1),gdc2(i,k,1)',gdc(i,k,1),gdc2(i,k,1) - endif ! !> - Calculate subsidence effect on clw ! - dsubclw=0. - dsubclwm=0. - dsubclws=0. +! dsubclw=0. +! dsubclwm=0. +! dsubclws=0. +! dp=100.*(p2d(i,k)-p2d(i,k+1)) +! if (clcw(i,k) .gt. -999.0 .and. clcw(i,k+1) .gt. -999.0 )then +! clwtot = cliw(i,k) + clcw(i,k) +! clwtot1= cliw(i,k+1) + clcw(i,k+1) +! dsubclw=((-edt(i)*zd(i,k+1)+zu(i,k+1))*clwtot1 & +! -(-edt(i)*zd(i,k) +zu(i,k)) *clwtot )*g/dp +! dsubclwm=((-edtm(i)*zdm(i,k+1)+zum(i,k+1))*clwtot1 & +! -(-edtm(i)*zdm(i,k) +zum(i,k)) *clwtot )*g/dp +! dsubclws=(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp +! dsubclw=dsubclw+(zu(i,k+1)*clwtot1-zu(i,k)*clwtot)*g/dp +! dsubclwm=dsubclwm+(zum(i,k+1)*clwtot1-zum(i,k)*clwtot)*g/dp +! dsubclws=dsubclws+(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp +! endif +! tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & +! +outqcm(i,k)*cutenm(i) & +! +dsubclw*xmb(i)+dsubclws*xmbs(i)+dsubclwm*xmbm(i) & +! ) +! tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) +! if (clcw(i,k) .gt. -999.0) then +! cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice +! clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water +! else +! cliw(i,k) = max(0.,cliw(i,k) + tem) +! endif +! +! enddo + +!> - FCT treats subsidence effect to cloud ice/water (begin) dp=100.*(p2d(i,k)-p2d(i,k+1)) + dtime_max=min(dtime_max,.5*dp) + po_cup(k)=.5*(p2d(i,k)+p2d(i,k+1)) if (clcw(i,k) .gt. -999.0 .and. clcw(i,k+1) .gt. -999.0 )then clwtot = cliw(i,k) + clcw(i,k) + if(clwtot.lt.1.e-32)clwtot=0. clwtot1= cliw(i,k+1) + clcw(i,k+1) - dsubclw=((-edt(i)*zd(i,k+1)+zu(i,k+1))*clwtot1 & - -(-edt(i)*zd(i,k) +zu(i,k)) *clwtot )*g/dp - dsubclwm=((-edtm(i)*zdm(i,k+1)+zum(i,k+1))*clwtot1 & - -(-edtm(i)*zdm(i,k) +zum(i,k)) *clwtot )*g/dp - dsubclws=(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp - dsubclw=dsubclw+(zu(i,k+1)*clwtot1-zu(i,k)*clwtot)*g/dp - dsubclwm=dsubclwm+(zum(i,k+1)*clwtot1-zum(i,k)*clwtot)*g/dp - dsubclws=dsubclws+(zus(i,k+1)*clwtot1-zus(i,k)*clwtot)*g/dp + if(clwtot1.lt.1.e-32)clwtot1=0. + clw_in1(k)=clwtot + massflx(k)=-(xmb(i) *( zu(i,k)- edt(i)* zd(i,k))) & + -(xmbm(i)*(zdm(i,k)-edtm(i)*zdm(i,k))) & + -(xmbs(i)*zus(i,k)) + trcflx_in1(k)=massflx(k)*.5*(clwtot+clwtot1) endif - tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & + enddo + + massflx (1)=0. + trcflx_in1(1)=0. + call fct1d3 (kstop,kte,dtime_max,po_cup, & + clw_in1,massflx,trcflx_in1,clw_ten1,g) + + do k=1,kstop + tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & +outqcm(i,k)*cutenm(i) & -! +dsubclw*xmb(i)+dsubclws*xmbs(i)+dsubclwm*xmbm(i) & - ) + +clw_ten1(k) & + ) tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) if (clcw(i,k) .gt. -999.0) then cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water - else + else cliw(i,k) = max(0.,cliw(i,k) + tem) - endif + endif - enddo - gdc(i,1,10)=forcing(i,1) - gdc(i,2,10)=forcing(i,2) - gdc(i,3,10)=forcing(i,3) - gdc(i,4,10)=forcing(i,4) - gdc(i,5,10)=forcing(i,5) - gdc(i,6,10)=forcing(i,6) - gdc(i,7,10)=forcing(i,7) - gdc(i,8,10)=forcing(i,8) - gdc(i,10,10)=xmb(i) - gdc(i,11,10)=xmbm(i) - gdc(i,12,10)=xmbs(i) - gdc(i,13,10)=hfx(i) - gdc(i,15,10)=qfx(i) - gdc(i,16,10)=pret(i)*3600. +! +!> calculate cloud water and cloud ice number concentrations +! + rho_dryar(i,k) = p2di(i,k)/(con_rd*t(i,k)) ! Density of dry air in kg m-3 + if (imp_physics == imp_physics_thompson) then + if ((tem*tem1)>1.e-5) then + gq0(i,k,ntinc) = max(0., gq0(i,k,ntinc) + & + make_IceNumber(tem*tem1*rho_dryar(i,k), t(i,k)) * & + (1/rho_dryar(i,k))) + end if + if ((tem*(1-tem1))>1.e-5) then + gq0(i,k,ntlnc) = max(0., gq0(i,k,ntlnc) + & + make_DropletNumber(tem*(1-tem1)*rho_dryar(i,k), nwfa(i,k)) & + * (1/rho_dryar(i,k))) + end if + end if + + enddo + + + gdc(i,1,10)=forcing(i,1) + gdc(i,2,10)=forcing(i,2) + gdc(i,3,10)=forcing(i,3) + gdc(i,4,10)=forcing(i,4) + gdc(i,5,10)=forcing(i,5) + gdc(i,6,10)=forcing(i,6) + gdc(i,7,10)=forcing(i,7) + gdc(i,8,10)=forcing(i,8) + gdc(i,10,10)=xmb(i) + gdc(i,11,10)=xmbm(i) + gdc(i,12,10)=xmbs(i) + gdc(i,13,10)=hfx(i) + gdc(i,15,10)=qfx(i) + gdc(i,16,10)=pret(i)*3600. if(ktop(i).gt.2 .and.pret(i).gt.0.)dt_mf(i,ktop(i)-1)=ud_mf(i,ktop(i)) endif enddo diff --git a/physics/cu_gf_driver.meta b/physics/cu_gf_driver.meta index 1969f9464..d3687a352 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -44,6 +44,14 @@ [ccpp-arg-table] name = cu_gf_driver_run type = scheme +[ntracer] + standard_name = number_of_tracers + long_name = number of tracers + units = count + dimensions = () + type = integer + intent = in + optional = F [garea] standard_name = cell_area long_name = grid cell area @@ -350,6 +358,65 @@ type = integer intent = in optional = F +[nwfa] + standard_name = water_friendly_aerosol_number_concentration + long_name = number concentration of water-friendly aerosols + units = kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[con_rd] + standard_name = gas_constant_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[gq0] + standard_name = tracer_concentration_updated_by_physics + long_name = tracer concentration updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension,number_of_tracers) + type = real + kind = kind_phys + intent = inout + optional = F +[ntinc] + standard_name = index_for_ice_cloud_number_concentration + long_name = tracer index for ice number concentration + units = index + dimensions = () + type = integer + intent = in + optional = F +[ntlnc] + standard_name = index_for_liquid_cloud_number_concentration + long_name = tracer index for liquid number concentration + units = index + dimensions = () + type = integer + intent = in + optional = F +[imp_physics] + standard_name = flag_for_microphysics_scheme + long_name = choice of microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imp_physics_thompson] + standard_name = flag_for_thompson_microphysics_scheme + long_name = choice of Thompson microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/docs/ccpp_doxyfile b/physics/docs/ccpp_doxyfile index 91c80c221..b435664e3 100644 --- a/physics/docs/ccpp_doxyfile +++ b/physics/docs/ccpp_doxyfile @@ -245,6 +245,10 @@ INPUT = pdftxt/mainpage.txt \ ../module_mp_thompson.F90 \ ../module_mp_radar.F90 \ ../mp_thompson_post.F90 \ +### HAFS + ../module_MP_FER_HIRES.F90 \ + ../mp_fer_hires.F90 \ + ../module_mp_fer_hires_pre.F90 \ ### utils ../funcphys.f90 \ ../physparam.f \ diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index c3da28334..269bf0b3a 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -596,6 +596,7 @@ subroutine drag_suite_run( & olss(i) = 0.0 ulow (i) = 0.0 dtfac(i) = 1.0 + rstoch(i) = 0.0 ldrag(i) = .false. icrilv(i) = .false. flag(i) = .true. diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 10533d99d..174e0c95c 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -26,7 +26,8 @@ end subroutine maximum_hourly_diagnostics_finalize !! #endif subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, & - imp_physics_gfdl, imp_physics_thompson, con_g, phil, & + imp_physics_gfdl, imp_physics_thompson, & + imp_physics_fer_hires,con_g, phil, & gt0, refl_10cm, refdmax, refdmax263k, u10m, v10m, & u10max, v10max, spd10max, pgr, t2m, q2m, t02max, & t02min, rh02max, rh02min, errmsg, errflg) @@ -34,7 +35,7 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, ! Interface variables integer, intent(in) :: im, levs logical, intent(in) :: reset, lradar - integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson + integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires real(kind_phys), intent(in ) :: con_g real(kind_phys), intent(in ) :: phil(im,levs) real(kind_phys), intent(in ) :: gt0(im,levs) @@ -66,9 +67,9 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, errflg = 0 !Calculate hourly max 1-km agl and -10C reflectivity - if (lradar .and. & - (imp_physics == imp_physics_gfdl .or. & - imp_physics == imp_physics_thompson)) then + if (lradar .and. (imp_physics == imp_physics_gfdl .or. & + imp_physics == imp_physics_thompson .or. & + imp_physics == imp_physics_fer_hires)) then allocate(refd(im)) allocate(refd263k(im)) call max_fields(phil,refl_10cm,con_g,im,levs,refd,gt0,refd263k) diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index df6f10913..5146ce2f0 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -57,6 +57,14 @@ type = integer intent = in optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [con_g] standard_name = gravitational_acceleration long_name = gravitational acceleration diff --git a/physics/module_MP_FER_HIRES.F90 b/physics/module_MP_FER_HIRES.F90 new file mode 100644 index 000000000..a736c640f --- /dev/null +++ b/physics/module_MP_FER_HIRES.F90 @@ -0,0 +1,2927 @@ +!>\file module_MP_FER_HIRES.F90 +!! "Modified" fer_hires microphysics - 11 July 2016 version +!! +! (1) Ice nucleation: Fletcher (1962) replaces Meyers et al. (1992) +! (2) Cloud ice is a simple function of the number concentration from (1), and it +! is no longer a fractional function of the large ice. Thus, the FLARGE & +! FSMALL parameters are no longer used. +! (3) T_ICE_init=-12 deg C provides a slight delay in the initial onset of ice. +! (4) NLImax is a function of rime factor (RF) and temperature. +! a) For RF>10, NLImax=1.e3. Mean ice diameters can exceed the 1 mm maximum +! size in the tables so that NLICE=NLImax=1.e3. +! b) Otherwise, NLImax is 10 L-1 at 0C and decreasing to 5 L-1 at <=-40C. +! NLICE>NLImax at the maximum ice diameter of 1 mm. +! (5) Can turn off ice processes by setting T_ICE & T_ICE_init to be < -100 deg C +! (6) Modified the homogeneous freezing of cloud water when TNLImax. +! (10) Ice deposition does not change the rime factor (RF) when RF>=10 & T>T_ICE. +! (11) Limit GAMMAS to <=1.5 (air resistance impact on ice fall speeds) +! (12) NSImax is maximum # conc of ice crsytals. At cold temperature NSImax is +! calculated based on assuming 10% of total ice content is due to cloud ice. +! +!-- Further modifications starting on 23 July 2015 +! (13) RHgrd is passed in as an input argument so that it can vary for different +! domains (RHgrd=0.98 for 12-km parent, 1.0 for 3-km nests) +! (14) Use the old "PRAUT" cloud water autoconversion *threshold* (QAUT0) + +!-- Further modifications starting on 28 July 2015 +! (15) Added calculations for radar reflectivity and number concentrations of +! rain (Nrain) and precipitating ice (Nsnow). +! (16) Removed double counting of air resistance term for riming onto ice (PIACW) +! (17) The maximum rime factor (RFmx) is now a function of MASSI(INDEXS), accounting +! for the increase in unrimed ice particle densities as values of INDEXS +! decrease from the maximum upper limit of 1000 microns to the lower limit of +! 50 microns, coinciding with the assumed size of cloud ice; see lines 1128-1134. +! (18) A new closure is used for updating the rime factor, which is described in +! detail near lines 1643-1682. The revised code is near lines 1683-1718. +! (19) Restructured the two-pass algorithm to be more robust, removed the HAIL +! & LARGE_RF logical variables so that NLICE>NLImax can occur. +! (20) Increased nsimax (see !aug27 below) +! (21) Modified the rain sedimentation (see two !aug27 blocks below) +! (22) NInuclei is the lower of Fletcher (1962), Cooper (1986), or NSImax. +! (23) NLImax is no longer used or enforced. Instead, INDEXS=MDImax when RF>20, +! else INDEXS is a function of temperature. Look for !sep10 comment. +! (24) An override was inserted for (18), such that the rime density is not diluted +! diluted when RF>20. Look for !sep10 comment. +! (25) Radar reflectivity calculations were changes to reduce radar bright bands, +! limit enhanced, mixed-phase reflectivity to RF>=20. Look for !sep10 comments. +! (26) NLICE is not to exceed NSI_max (250 L^-1) when RF<20. Look for !sep16 comments. +! Commented out! (28) Increase hail fall speeds using Thompson et al. (2008). Look for !sep22 comments. +! (29) Modify NLImax, INDEXS for RF>=20. Look for !sep22 comments. +! (30) Check on NSmICE, Vci based on whether FLIMASS<1. Look for !sep22a comments. +! Revised in (34)! (31) Introduced RFlag logical, which if =T enforces a lower limit of drop sizes not +! to go below INDEXRmin and N0r is adjusted. Look for !nov25 comments (corrections, +! refinements to sep25 & nov18 versions, includes an additional fix in nov25-fix). +! Also set INDEXRmin=500 rather than 250 microns. +!----------------------------------------------------------------------------- +!--- The following changes now refer to dates when those were made in 2016. +!----------------------------------------------------------------------------- +! (32) Convective (RF>=20, Ng~10 L^-1, RHOg~500 kg m^-3), transition (RF=10, Ng~25 L^-1, +! RHOg~300 kg m^-3), & stratiform (RF<2) profiles are blended based on RF. !mar08 +! (33) Fixed bug in Biggs' freezing, put back in collisional drop freezing. !mar03 +! (34) Changes in (31) are revised so that INDEXRmin at and below 0C level is +! based on a rain rate equal to the snowfall rate above the 0C level. !mar03 +! (35) Increase radar reflectivity when RF>10 and RQSnew > 2.5 g m^-3. !mar12 +! (36) !mar10 combines all elements of (32)-(35) together. +! (37) Bug fixes for the changes in (34) and the RFLAG variable !apr18 +! (38) Revised Schumann-Ludlam limit. !apr18 +! (39) Simplified PCOND (cloud cond/evap) calculation !apr21 +! (40) Slight change in calculating RF. !apr22 +! (41) Reduce RF values for calculating mean sizes of snow, graupel, sleet/hail !apr22a +! (42) Increase reflectivity from large, wet, high rime factor ice (graupel) by +! assuming |Kw|**2/|Ki|**2 = 0.224 (Smith, 1984, JCAM). +! (43) Major restructuring of code to allow N0r to vary from N0r0 !may11 +! (44) More major restructuring of code to use fixed XLS, XLV, XLF !may12 +! (45) Increased VEL_INC ~ VrimeF**2, put the enhanced graupel/hail fall speeds +! from Thompson into the code but only in limited circumstances, restructured +! and streamlined the INDEXS calculation, removed the upper limit for +! for the vapor mixing ratio is at water saturation when calculating ice +! deposition, and N0r is gradually increased for conditions supporting +! drizzle when rain contents decrease below 0.25 g/m**3. !may17 +! (46) The may11 code changes that increase N0r0 when rain contents exceed 1 g m^-3 +! have been removed, limit the number of iterations calculating final rain +! parameters, remove the revised N0r calculation for reflectivity. All of +! the changes following those made in the may10 code. !may20 +! (47) Reduce the assumed # concentration of hail/sleet when RF>10 from 5 L^-1 to +! 1 L^-1, and also reduce it for graupel when RF>5 from 10 L^-1 to 5 L^-1. +! This is being done to try and make greater use of the Thompson graupel/hail +! fallspeeds by having INDEXS==MDImax. +! (48) Increased NCW from 200e6 to 300e6 for a more delayed onset of drizzle, +! simplified drizzle algorithm to reduce/eliminate N0r bulls eyes and to allow +! for supercooled drizzle, and set limits for 8.e6 <= N0r (m^-4) <= 1.e9 !may31 +! (49) Further restructuring of code to better define STRAT, DRZL logicals, +! add these rain flags to mprates arrays !jun01 +! (50) Increase in reflectivity due to wet ice was commented out. +! (51) Fixed minor bug to update INDEXR2 in the "rain_pass: do" loop. !jun13 +! (52) Final changes to Nsnow for boosting reflectivities from ice for +! mass contents exceeding 5 g m^-3. !jun16 +! (53) Cosmetic changes only that do not affect the calculations. Removed old, unused +! diagnostic arrays. Updated comments. +! +!----------------------------------------------------------------------------- +! + MODULE MODULE_MP_FER_HIRES +! +!----------------------------------------------------------------------------- + +#ifdef MPI + USE mpi +#endif + USE machine +!MZ +!MZ USE MODULE_CONSTANTS,ONLY : PI, CP, EPSQ, GRAV=>G, RHOL=>RHOWATER, & +!MZ RD=>R_D, RV=>R_V, T0C=>TIW, EPS=>EP_2, EPS1=>EP_1, CLIQ, CICE, & +!MZ XLV +!MZ +!MZ temporary values copied from module_CONSTANTS; ideally they come from host model +!side + REAL, PARAMETER :: pi=3.141592653589793 ! ludolf number + REAL, PARAMETER :: cp=1004.6 ! spec. heat for dry air at constant pressure + REAL, PARAMETER :: epsq=1.e-12 ! floor value for specific humidity (kg/kg) + REAL, PARAMETER :: grav= 9.8060226 ! gravity + REAL, PARAMETER :: RHOL=1000. ! density of water (kg/m3) + REAL, PARAMETER :: RD=287.04 ! gas constant for dry air + REAL, PARAMETER :: RV=461.6 ! gas constant for water vapor + REAL, PARAMETER :: T0C= 273.15 ! melting point + REAL, PARAMETER :: EPS=RD/RV + REAL, PARAMETER :: EPS1=RV/RD-1. + REAL, PARAMETER :: CLIQ = 4190. ! MZ: inconsistent value below + REAL, PARAMETER :: CICE = 2106. + REAL, PARAMETER :: XLV = 2.5E6 +!----------------------------------------------------------------------------- + PUBLIC :: FERRIER_INIT_HR, GPVS_HR,FPVS,FPVS0,NX +!----------------------------------------------------------------------------- + REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, C_NR, Crain, & !jul28 + & CRACW, ARAUT, BRAUT, ESW0, RFmx1, ARcw, RH_NgC, RH_NgT, & !jul31 !mar08 + & RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DR4, RR_DR5, RR_DRmax, & !may17 + & BETA6, & + & RQhail, AVhail, BVhail, QAUT0 !may17 +! + INTEGER,PRIVATE,PARAMETER :: INDEXRstrmax=500 !mar03, stratiform maximum + REAL,PUBLIC,SAVE :: CN0r0, CN0r_DMRmin, CN0r_DMRmax, & + RFmax, RQR_DRmax, RQR_DRmin +! + INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35 + REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH_NMM +! + REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, & + & DelDMI=1.e-6,XMImin=1.e6*DMImin + REAL, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536 + INTEGER, PUBLIC,PARAMETER :: MDImin=XMImin, MDImax=XMImax + REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: & + & ACCRI,VSNOWI,VENTI1,VENTI2 + REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: SDENS !-- For RRTM +! + REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=1.0e-3, & + & DelDMR=1.e-6, XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax + INTEGER, PUBLIC,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax +! + REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: & + & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2 +! + INTEGER, PRIVATE,PARAMETER :: Nrime=40 + REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF +! + INTEGER,PARAMETER :: NX=7501 + REAL, PARAMETER :: XMIN=180.0,XMAX=330.0 + REAL, DIMENSION(NX),PUBLIC,SAVE :: TBPVS,TBPVS0 + REAL, PUBLIC,SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS +! + REAL,DIMENSION(MY_T2+8) :: MP_RESTART_STATE + REAL,DIMENSION(nx) :: TBPVS_STATE,TBPVS0_STATE +! + REAL, PRIVATE,PARAMETER :: CVAP=1846., XLF=3.3358e+5, XLS=XLV+XLF & + & ,EPSQ1=1.001*EPSQ, RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV & + & ,RRHOL=1./RHOL, XLV1=XLV/CP, XLF1=XLF/CP, XLS1=XLS/CP & + & ,XLV2=XLV*XLV/RV, XLS2=XLS*XLS/RV & + & ,XLV3=XLV*XLV*RCPRV, XLS3=XLS*XLS*RCPRV & +!--- Constants specific to the parameterization follow: +!--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation + & ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT & + & ,C1=1./3. & + & ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, DMR4=0.45E-3 & + & ,DMR5=0.67E-3 & + & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3 & + & ,XMR4=1.e6*DMR4, XMR5=1.e6*DMR5, RQRmix=0.05E-3, RQSmix=1.E-3 & !jul28 !apr27 + & ,Cdry=1.634e13, Cwet=1./.224 !jul28 !apr27 + INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3, MDR4=XMR4 & + & , MDR5=XMR5 + +!-- Debug 20120111 +LOGICAL, SAVE :: WARN1=.TRUE.,WARN2=.TRUE.,WARN3=.TRUE.,WARN5=.TRUE. +REAL, SAVE :: Pwarn=75.E2, QTwarn=1.E-3 +INTEGER, PARAMETER :: MAX_ITERATIONS=10 + +! +! ====================================================================== +!--- Important tunable parameters that are exported to other modules +! * T_ICE - temperature (C) threshold at which all remaining liquid water +! is glaciated to ice +! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs +! +!-- To turn off ice processes, set T_ICE & T_ICE_init to <= -100. (i.e., -100 C) +! +! * NSImax - maximum number concentrations (m**-3) of small ice crystals +! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet) +! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 1.0 mm +! * N0rmin - minimum intercept (m**-4) for rain drops +! * NCW - number concentrations of cloud droplets (m**-3) +! ====================================================================== + REAL, PUBLIC,PARAMETER :: & + & RHgrd_in=1. & + &, P_RHgrd_out=850.E2 & + & ,T_ICE=-40. & + & ,T_ICEK=T0C+T_ICE & + & ,T_ICE_init=-12. & + & ,NSI_max=250.E3 & + & ,NLImin=1.0E3 & + & ,N0r0=8.E6 & + & ,N0rmin=1.E4 & +!! based on Aligo's email,NCW is changed to 250E6 + & ,NCW=250.E6 + !HWRF & ,NCW=300.E6 !- 100.e6 (maritime), 500.e6 (continental) + +!--- Other public variables passed to other routines: + REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI +! + + CONTAINS +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& +!----------------------------------------------------------------------- + +!>\ingroup hafs_famp +!! This is the driver scheme of Ferrier-Aligo microphysics scheme. +!! NOTE: The only differences between FER_HIRES and FER_HIRES_ADVECT +!! is that the QT, and F_* are all local variables in the advected +!! version, and QRIMEF is only in the advected version. The innards +!! are all the same. + SUBROUTINE FER_HIRES (DT,RHgrd, & + & dz8w,rho_phy,p_phy,pi_phy,th_phy,t_phy, & + & q,qt, & + & LOWLYR,SR, & + & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, & + & QC,QR,QS, & + & RAINNC,RAINNCV, & + & threads, & + & ims,ime, jms,jme, lm, & + & d_ss, & + & refl_10cm,DX1 ) +!----------------------------------------------------------------------- + IMPLICIT NONE +!----------------------------------------------------------------------- + INTEGER,INTENT(IN) :: D_SS,IMS,IME,JMS,JME,LM,DX1 + REAL, INTENT(IN) :: DT,RHgrd + INTEGER, INTENT(IN) :: THREADS + REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme, lm):: & + & dz8w,p_phy,pi_phy,rho_phy + REAL, INTENT(INOUT), DIMENSION(ims:ime, jms:jme, lm):: & + & th_phy,t_phy,q,qt + REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme, lm ) :: & + & qc,qr,qs + REAL, INTENT(INOUT), DIMENSION(ims:ime, jms:jme,lm) :: & + & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY + REAL, INTENT(OUT), DIMENSION(ims:ime, jms:jme,lm) :: & + & refl_10cm + REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: & + & RAINNC,RAINNCV + REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR +! + INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR + +!----------------------------------------------------------------------- +! LOCAL VARS +!----------------------------------------------------------------------- + +! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related +! the microphysics scheme. Instead, they will be used by Eta precip +! assimilation. + + REAL, DIMENSION( ims:ime, jms:jme,lm ) :: & + & TLATGS_PHY,TRAIN_PHY + REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC + + INTEGER :: I,J,K,KK + REAL :: wc +!------------------------------------------------------------------------ +! For subroutine EGCP01COLUMN_hr +!----------------------------------------------------------------------- + INTEGER :: LSFC,I_index,J_index,L + INTEGER,DIMENSION(ims:ime,jms:jme) :: LMH + REAL :: TC,QI,QRdum,QW,Fice,Frain,DUM,ASNOW,ARAIN + REAL,DIMENSION(lm) :: P_col,Q_col,T_col,WC_col, & + RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL,pcond1d, & + pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d, & + pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col, & + NR_col,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d, & + INDEXS1d,INDEXR1d,RFlag1d,RHC_col +! +!----------------------------------------------------------------------- +!********************************************************************** +!----------------------------------------------------------------------- +! + +! MZ: HWRF practice start +!---------- +!2015-03-30, recalculate some constants which may depend on phy time step + CALL MY_GROWTH_RATES_NMM_hr (DT) + +!--- CIACW is used in calculating riming rates +! The assumed effective collection efficiency of cloud water rimed onto +! ice is =0.5 below: +! + CIACW=DT*0.25*PI*0.5*(1.E5)**C1 +! +!--- CIACR is used in calculating freezing of rain colliding with large ice +! The assumed collection efficiency is 1.0 +! + CIACR=PI*DT +! +!--- CRACW is used in calculating collection of cloud water by rain (an +! assumed collection efficiency of 1.0) +! + CRACW=DT*0.25*PI*1.0 +! +!-- See comments in subroutine etanewhr_init starting with variable RDIS= +! + BRAUT=DT*1.1E10*BETA6/NCW + + !write(*,*)'dt=',dt + !write(*,*)'pi=',pi + !write(*,*)'c1=',c1 + !write(*,*)'ciacw=',ciacw + !write(*,*)'ciacr=',ciacr + !write(*,*)'cracw=',cracw + !write(*,*)'araut=',araut + !write(*,*)'braut=',braut +!! END OF adding, 2015-03-30 +!----------- +! MZ: HWRF practice end +! + + DO j = jms,jme + DO i = ims,ime + ACPREC(i,j)=0. + APREC (i,j)=0. + PREC (i,j)=0. + SR (i,j)=0. + ENDDO + DO k = 1,lm + DO i = ims,ime + TLATGS_PHY (i,j,k)=0. + TRAIN_PHY (i,j,k)=0. + ENDDO + ENDDO + ENDDO + +!----------------------------------------------------------------------- +!-- Start of original driver for EGCP01COLUMN_hr +!----------------------------------------------------------------------- +! + DO J=JMS,JME + DO I=IMS,IME + LSFC=LM-LOWLYR(I,J)+1 ! "L" of surface + DO K=1,LM + DPCOL(K)=RHO_PHY(I,J,K)*GRAV*dz8w(I,J,K) + ENDDO +! +!--- Initialize column data (1D arrays) +! + L=LM +!-- qt = CWM, total condensate + IF (qt(I,J,L) .LE. EPSQ) qt(I,J,L)=EPSQ + F_ice_phy(I,J,L)=1. + F_rain_phy(I,J,L)=0. + F_RimeF_phy(I,J,L)=1. + do L=LM,1,-1 +! +!--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop +! + P_col(L)=P_phy(I,J,L) +! +!--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc) +! + THICK_col(L)=DPCOL(L)*RGRAV + T_col(L)=T_phy(I,J,L) + TC=T_col(L)-T0C + Q_col(L)=max(EPSQ, q(I,J,L)) + IF (qt(I,J,L) .LE. EPSQ1) THEN + WC_col(L)=0. + IF (TC .LT. T_ICE) THEN + F_ice_phy(I,J,L)=1. + ELSE + F_ice_phy(I,J,L)=0. + ENDIF + F_rain_phy(I,J,L)=0. + F_RimeF_phy(I,J,L)=1. + ELSE + WC_col(L)=qt(I,J,L) + +!-- Debug 20120111 +! TC==TC will fail if NaN, preventing unnecessary error messages +IF (WC_col(L)>QTwarn .AND. P_col(L)1 g/kg condensate in stratosphere; I,J,L,TC,P,QT=', & + I,J,L,TC,.01*P_col(L),1000.*WC_col(L) + QTwarn=MAX(WC_col(L),10.*QTwarn) + Pwarn=MIN(P_col(L),0.5*Pwarn) +ENDIF +!-- TC/=TC will pass if TC is NaN +IF (WARN5 .AND. TC/=TC) THEN + WRITE(0,*) 'WARN5: NaN temperature; I,J,L,P=',I,J,L,.01*P_col(L) + WARN5=.FALSE. +ENDIF + + ENDIF + IF (T_ICE<=-100.) F_ice_phy(I,J,L)=0. +! ! +! !--- Determine composition of condensate in terms of +! ! cloud water, ice, & rain +! ! + WC=WC_col(L) + QI=0. + QRdum=0. + QW=0. + Fice=F_ice_phy(I,J,L) + Frain=F_rain_phy(I,J,L) +! + IF (Fice .GE. 1.) THEN + QI=WC + ELSE IF (Fice .LE. 0.) THEN + QW=WC + ELSE + QI=Fice*WC + QW=WC-QI + ENDIF +! + IF (QW.GT.0. .AND. Frain.GT.0.) THEN + IF (Frain .GE. 1.) THEN + QRdum=QW + QW=0. + ELSE + QRdum=Frain*QW + QW=QW-QRdum + ENDIF + ENDIF + IF (QI .LE. 0.) F_RimeF_phy(I,J,L)=1. + RimeF_col(L)=F_RimeF_phy(I,J,L) ! (real) + QI_col(L)=QI + QR_col(L)=QRdum + QW_col(L)=QW +!GFDL => New. Added RHC_col to allow for height- and grid-dependent values for +!GFDL the relative humidity threshold for condensation ("RHgrd") +!6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa +!------------------------------------------------------------ + IF(DX1 .GE. 10 .AND. P_col(L)0) associated with snow +! + APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying + PREC(I,J)=PREC(I,J)+APREC(I,J) + ACPREC(I,J)=ACPREC(I,J)+APREC(I,J) + IF(APREC(I,J) .LT. 1.E-8) THEN + SR(I,J)=0. + ELSE + SR(I,J)=RRHOL*ASNOW/APREC(I,J) + ENDIF +! +!####################################################################### +!####################################################################### +! + enddo ! End "I" loop + enddo ! End "J" loop +! +!----------------------------------------------------------------------- +!-- End of original driver for EGCP01COLUMN_hr +!----------------------------------------------------------------------- +! + DO j = jms,jme + do k = lm, 1, -1 + DO i = ims,ime + th_phy(i,j,k) = t_phy(i,j,k)/pi_phy(i,j,k) + WC=qt(I,J,K) + QS(I,J,K)=0. + QR(I,J,K)=0. + QC(I,J,K)=0. +! + IF(F_ICE_PHY(I,J,K)>=1.)THEN + QS(I,J,K)=WC + ELSEIF(F_ICE_PHY(I,J,K)<=0.)THEN + QC(I,J,K)=WC + ELSE + QS(I,J,K)=F_ICE_PHY(I,J,K)*WC + QC(I,J,K)=WC-QS(I,J,K) + ENDIF +! + IF(QC(I,J,K)>0..AND.F_RAIN_PHY(I,J,K)>0.)THEN + IF(F_RAIN_PHY(I,J,K).GE.1.)THEN + QR(I,J,K)=QC(I,J,K) + QC(I,J,K)=0. + ELSE + QR(I,J,K)=F_RAIN_PHY(I,J,K)*QC(I,J,K) + QC(I,J,K)=QC(I,J,K)-QR(I,J,K) + ENDIF + ENDIF + ENDDO !- i + ENDDO !- k + ENDDO !- j +! +!- Update rain (convert from m to kg/m**2, which is also equivalent to mm depth) +! + DO j=jms,jme + DO i=ims,ime + RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j) + RAINNCV(i,j)=APREC(i,j)*1000. + ENDDO + ENDDO +! +!----------------------------------------------------------------------- +! + END SUBROUTINE FER_HIRES +! +!----------------------------------------------------------------------- +! +!############################################################################### +! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL +! (1) Represents sedimentation by preserving a portion of the precipitation +! through top-down integration from cloud-top. Modified procedure to +! Zhao and Carr (1997). +! (2) Microphysical equations are modified to be less sensitive to time +! steps by use of Clausius-Clapeyron equation to account for changes in +! saturation mixing ratios in response to latent heating/cooling. +! (3) Prevent spurious temperature oscillations across 0C due to +! microphysics. +! (4) Uses lookup tables for: calculating two different ventilation +! coefficients in condensation and deposition processes; accretion of +! cloud water by precipitation; precipitation mass; precipitation rate +! (and mass-weighted precipitation fall speeds). +! (5) Assumes temperature-dependent variation in mean diameter of large ice +! (Houze et al., 1979; Ryan et al., 1996). +! -> 8/22/01: This relationship has been extended to colder temperatures +! to parameterize smaller large-ice particles down to mean sizes of MDImin, +! which is 50 microns reached at -55.9C. +! (6) Attempts to differentiate growth of large and small ice, mainly for +! improved transition from thin cirrus to thick, precipitating ice +! anvils. +! (7) Top-down integration also attempts to treat mixed-phase processes, +! allowing a mixture of ice and water. Based on numerous observational +! studies, ice growth is based on nucleation at cloud top & +! subsequent growth by vapor deposition and riming as the ice particles +! fall through the cloud. There are two modes of ice nucleation +! following Meyers et al. (JAM, 1992): +! a) Deposition & condensation freezing nucleation - eq. (2.4) when +! air is supersaturated w/r/t ice +! b) Contact freezing nucleation - eq. (2.6) in presence of cloud water +! (8) Depositional growth of newly nucleated ice is calculated for large time +! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals +! using their ice crystal masses calculated after 600 s of growth in water +! saturated conditions. The growth rates are normalized by time step +! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993). +! (9) Ice precipitation rates can increase due to increase in response to +! cloud water riming due to (a) increased density & mass of the rimed +! ice, and (b) increased fall speeds of rimed ice. +!############################################################################### +!############################################################################### +! +!>\ingroup hafs_famp +!! This is the grid-scale microphysical processes of Ferrier-Aligo microphysics +!! scheme (i.e., condensation and precipitation). +!!\param arain accumulated rainfall at the surface (kg) +!!\param asnow accumulated snowfall at the surface (kg) +!!\param dtph physics time step (s) +!!\param rhc_col vertical column of threshold relative humidity for onset of +!! condensation (ratio) +!!\param i_index i index +!!\param j_index j index +!!\param lsfc Eta level of level above surface, ground +!!\param p_col vertical column of model pressure (Pa) +!!\param qi_col vertical column of model ice mixing ratio (kg/kg) +!!\param qr_col vertical column of model rain ratio (kg/kg) +!!\param q_col vertical column of model water vapor specific humidity (kg/kg) +!!\param qw_col +!!\param rimef_col +!!\param t_col +!!\param thick_col +!!\param wc_col +!!\param lm +!!\param pcond1d +!!\param pidep1d +!!\param piacw1d +!!\param piacwi1d + SUBROUTINE EGCP01COLUMN_hr ( ARAIN, ASNOW, DTPH, RHC_col, & + & I_index, J_index, LSFC, & + & P_col, QI_col, QR_col, Q_col, QW_col, RimeF_col, T_col, & + & THICK_col, WC_col ,LM,pcond1d,pidep1d, & + & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, & + & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, & + & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28 + & RFlag1d,DX1) !jun01 +! +!############################################################################### +!############################################################################### +! +!------------------------------------------------------------------------------- +!----- NOTE: Code is currently set up w/o threading! +!------------------------------------------------------------------------------- +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation +! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001 +! PRGRMMR: Jin (Modification for WRF structure) +!------------------------------------------------------------------------------- +! ABSTRACT: +! * Merges original GSCOND & PRECPD subroutines. +! * Code has been substantially streamlined and restructured. +! * Exchange between water vapor & small cloud condensate is calculated using +! the original Asai (1965, J. Japan) algorithm. See also references to +! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al. +! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR) +! parameterization. +!------------------------------------------------------------------------------- +! +! USAGE: +! * CALL EGCP01COLUMN_hr FROM SUBROUTINE EGCP01DRV +! +! INPUT ARGUMENT LIST: +! DTPH - physics time step (s) +! RHgrd - threshold relative humidity (ratio) for onset of condensation +! I_index - I index +! J_index - J index +! LSFC - Eta level of level above surface, ground +! P_col - vertical column of model pressure (Pa) +! QI_col - vertical column of model ice mixing ratio (kg/kg) +! QR_col - vertical column of model rain ratio (kg/kg) +! Q_col - vertical column of model water vapor specific humidity (kg/kg) +! QW_col - vertical column of model cloud water mixing ratio (kg/kg) +! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below) +! T_col - vertical column of model temperature (deg K) +! THICK_col - vertical column of model mass thickness (density*height increment) +! WC_col - vertical column of model mixing ratio of total condensate (kg/kg) +! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL +! +! +! OUTPUT ARGUMENT LIST: +! ARAIN - accumulated rainfall at the surface (kg) +! ASNOW - accumulated snowfall at the surface (kg) +! Q_col - vertical column of model water vapor specific humidity (kg/kg) +! WC_col - vertical column of model mixing ratio of total condensate (kg/kg) +! QW_col - vertical column of model cloud water mixing ratio (kg/kg) +! QI_col - vertical column of model ice mixing ratio (kg/kg) +! QR_col - vertical column of model rain ratio (kg/kg) +! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below) +! T_col - vertical column of model temperature (deg K) +! DBZ_col - vertical column of radar reflectivity (dBZ) +! NR_col - vertical column of rain number concentration (m^-3) +! NS_col - vertical column of snow number concentration (m^-3) +! +! OUTPUT FILES: +! NONE +! +! Subprograms & Functions called: +! * Real Function CONDENSE - cloud water condensation +! * Real Function DEPOSIT - ice deposition (not sublimation) +! * Integer Function GET_INDEXR - estimate the mean size of raindrops (microns) +! +! UNIQUE: NONE +! +! LIBRARY: NONE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +! MACHINE : IBM SP +! +!------------------------------------------------------------------------- +!--------------- Arrays & constants in argument list --------------------- +!------------------------------------------------------------------------- +! + IMPLICIT NONE +! + INTEGER,INTENT(IN) :: LM,I_index, J_index, LSFC,DX1 + REAL,INTENT(IN) :: DTPH + REAL,INTENT(INOUT) :: ARAIN, ASNOW + REAL,DIMENSION(LM),INTENT(INOUT) :: P_col, QI_col,QR_col & + & ,Q_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col,pcond1d & + & ,pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d & + & ,pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col,NR_col & + & ,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d & !jun01 + & ,INDEXR1d,RFlag1d,RHC_col !jun01 +! +!-------------------------------------------------------------------------------- +!--- The following arrays are integral calculations based on the mean +! snow/graupel diameters, which vary from 50 microns to 1000 microns +! (1 mm) at 1-micron intervals and assume exponential size distributions. +! The values are normalized and require being multipled by the number +! concentration of large ice (NLICE). +!--------------------------------------- +! - VENTI1 - integrated quantity associated w/ ventilation effects +! (capacitance only) for calculating vapor deposition onto ice +! - VENTI2 - integrated quantity associated w/ ventilation effects +! (with fall speed) for calculating vapor deposition onto ice +! - ACCRI - integrated quantity associated w/ cloud water collection by ice +! - MASSI - integrated quantity associated w/ ice mass +! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate +! precipitation rates +! - VEL_RF - velocity increase of rimed particles as functions of crude +! particle size categories (at 0.1 mm intervals of mean ice particle +! sizes) and rime factor (different values of Rime Factor of 1.1**N, +! where N=0 to Nrime). +!-------------------------------------------------------------------------------- +!--- The following arrays are integral calculations based on the mean +! rain diameters, which vary from 50 microns to 1000 microns +! (1 mm) at 1-micron intervals and assume exponential size distributions. +! The values are normalized and require being multiplied by the rain intercept +! (N0r). +!--------------------------------------- +! - VENTR1 - integrated quantity associated w/ ventilation effects +! (capacitance only) for calculating evaporation from rain +! - VENTR2 - integrated quantity associated w/ ventilation effects +! (with fall speed) for calculating evaporation from rain +! - ACCRR - integrated quantity associated w/ cloud water collection by rain +! - MASSR - integrated quantity associated w/ rain +! - VRAIN - mass-weighted fall speed of rain, used to calculate +! precipitation rates +! - RRATE - precipitation rates, which should also be equal to RHO*QR*VRAIN +! +!------------------------------------------------------------------------- +!------- Key parameters, local variables, & important comments --------- +!----------------------------------------------------------------------- +! +!--- TOLER => Tolerance or precision for accumulated precipitation +! + REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, & + Xratio=.025, Zmin=0.01, DBZmin=-20. +! +!--- If BLEND=1: +! precipitation (large) ice amounts are estimated at each level as a +! blend of ice falling from the grid point above and the precip ice +! present at the start of the time step (see TOT_ICE below). +!--- If BLEND=0: +! precipitation (large) ice amounts are estimated to be the precip +! ice present at the start of the time step. +! +!--- Extended to include sedimentation of rain on 2/5/01 +! + REAL, PARAMETER :: BLEND=1. +! +!--- This variable is for debugging purposes (if .true.) +! + LOGICAL, PARAMETER :: PRINT_diag=.false. +! +!----------------------------------------------------------------------- +!--- Local variables +!----------------------------------------------------------------------- +! + REAL :: EMAIRI, N0r, NLICE, NSmICE, NInuclei, Nrain, Nsnow, Nmix + REAL :: RHgrd + LOGICAL :: CLEAR, ICE_logical, DBG_logical, RAIN_logical, & + STRAT, DRZL + INTEGER :: INDEX_MY,INDEXR,INDEXR1,INDEXR2,INDEXS,IPASS,ITDX,IXRF,& + & IXS,LBEF,L,INDEXRmin,INDEXS0C,IDR !mar03 !may20 +! +! + REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, & + & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, & + & DENOMI,DENOMW,DENOMWI,DIDEP, & + & DIEVP,DIFFUS,DLI,DTRHO,DUM,DUM1,DUM2,DUM3, & + & DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLIMASS, & + & FWR,FWS,GAMMAR,GAMMAS, & + & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, & + & PIEVP,PILOSS,PIMLT,PINIT,PP,PRACW,PRAUT,PREVP,PRLOSS, & + & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, & + & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,Q,QW,QWnew,Rcw, & + & RFACTOR,RFmx,RFrime,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, & + & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, & + & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, & + & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, & + & VSNOW1,WC,WCnew,WSgrd,WS,WSnew,WV,WVnew, & + & XLI,XLIMASS,XRF, & + & NSImax,QRdum,QSmICE,QLgIce,RQLICE,VCI,TIMLT, & + & RQSnew,RQRnew,Zrain,Zsnow,Ztot,RHOX0C,RFnew,PSDEP,DELS !mar03 !apr22 + REAL, SAVE :: Revised_LICE=1.e-3 !-- kg/m**3 +! +!####################################################################### +!########################## Begin Execution ############################ +!####################################################################### +! +! + ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2) + VRAIN1=0. ! Rain fall speeds into grib box from above (m/s) + VSNOW1=0. ! Ice fall speeds into grib box from above (m/s) + ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2) + INDEXS0C=MDImin ! Mean snow/graupel diameter just above (<0C) freezing level (height) + RHOX0C=22.5 ! Estimated ice density at 0C (kg m^-3) !mar03 + TIMLT=0. ! Total ice melting in a layer (drizzle detection) + STRAT=.FALSE. ! Stratiform rain DSD below melting level !may11 + DRZL=.FALSE. ! Drizzle DSD below melting level !may23 +! +!----------------------------------------------------------------------- +!------------ Loop from top (L=1) to surface (L=LSFC) ------------------ +!----------------------------------------------------------------------- +! +big_loop: DO L=LM,1,-1 + pcond1d(L)=0. + pidep1d(L)=0. + piacw1d(L)=0. + piacwi1d(L)=0. + piacwr1d(L)=0. + piacr1d(L)=0. + picnd1d(L)=0. + pievp1d(L)=0. + pimlt1d(L)=0. + praut1d(L)=0. + pracw1d(L)=0. + prevp1d(L)=0. + pisub1d(L)=0. + pevap1d(L)=0. + vsnow1d(L)=0. + vrain11d(L)=0. + vrain21d(L)=0. + vci1d(L)=0. + NSmICE1d(L)=0. + DBZ_col(L)=DBZmin + NR_col(L)=0. + NS_col(L)=0. + INDEXR1d(L)=0. + INDEXS1d(L)=0. + RFlag1d(L)=0. !jun01 +! +!--- Skip this level and go to the next lower level if no condensate +! and very low specific humidities +! +!--- Check if any rain is falling into layer from above +! + IF (ARAIN .GT. CLIMIT) THEN + CLEAR=.FALSE. + VRAIN1=0. + ELSE + CLEAR=.TRUE. + ARAIN=0. + ENDIF +! +!--- Check if any ice is falling into layer from above +! +!--- NOTE that "SNOW" in variable names is often synonomous with +! large, precipitation ice particles +! + IF (ASNOW .GT. CLIMIT) THEN + CLEAR=.FALSE. + VSNOW1=0. + ELSE + ASNOW=0. + ENDIF +! +!----------------------------------------------------------------------- +!------------ Proceed with cloud microphysics calculations ------------- +!----------------------------------------------------------------------- +! + TK=T_col(L) ! Temperature (deg K) + TC=TK-T0C ! Temperature (deg C) + PP=P_col(L) ! Pressure (Pa) + Q=Q_col(L) ! Specific humidity of water vapor (kg/kg) + WV=Q/(1.-Q) ! Water vapor mixing ratio (kg/kg) + WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg) + RHgrd=RHC_col(L) ! Threshold relative humidity for the onset of condensation +! +!----------------------------------------------------------------------- +!--- Moisture variables below are mixing ratios & not specifc humidities +!----------------------------------------------------------------------- +! +!--- This check is to determine grid-scale saturation when no condensate is present +! + ESW=MIN(1000.*FPVS0(TK),0.99*PP) ! Saturation vapor pressure w/r/t water + QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water + WS=QSW ! General saturation mixing ratio (water/ice) + QSI=QSW ! Saturation mixing ratio w/r/t ice + IF (TC .LT. 0.) THEN + ESI=MIN(1000.*FPVS(TK),0.99*PP) ! Saturation vapor pressure w/r/t ice + QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water + WS=QSI ! General saturation mixing ratio (water/ice) + ENDIF +! +!--- Effective grid-scale Saturation mixing ratios +! + QSWgrd=RHgrd*QSW + QSIgrd=RHgrd*QSI + WSgrd=RHgrd*WS +! +!--- Check if air is subsaturated and w/o condensate +! + IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE. +! +!----------------------------------------------------------------------- +!-- Loop to the end if in clear, subsaturated air free of condensate --- +!----------------------------------------------------------------------- +! + IF (CLEAR) THEN + STRAT=.FALSE. !- Reset stratiform rain flag + DRZL=.FALSE. !- Reset drizzle flag + INDEXRmin=MDRmin !- Reset INDEXRmin + TIMLT=0. !- Reset accumulated ice melting + CYCLE big_loop + ENDIF +! +!----------------------------------------------------------------------- +!--------- Initialize RHO, THICK & microphysical processes ------------- +!----------------------------------------------------------------------- +! +! +!--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity; +! (see pp. 63-65 in Fleagle & Businger, 1963) +! + RHO=PP/(RD*TK*(1.+EPS1*Q)) ! Air density (kg/m**3) + RRHO=1./RHO ! Reciprocal of air density + DTRHO=DTPH*RHO ! Time step * air density + BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density + THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc) +! + ARAINnew=0. ! Updated accumulated rainfall + ASNOWnew=0. ! Updated accumulated snowfall + QI=QI_col(L) ! Ice mixing ratio + QInew=0. ! Updated ice mixing ratio + QR=QR_col(L) ! Rain mixing ratio + QRnew=0. ! Updated rain ratio + QW=QW_col(L) ! Cloud water mixing ratio + QWnew=0. ! Updated cloud water ratio +! + PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg) + PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg) + PINIT=0. ! Ice initiation (part of PIDEP calculation, kg/kg) + PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0) + PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0) + PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0) + PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0) + PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg) + PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg) + PIMLT=0. ! Melting ice (kg/kg; >0) + PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0) + PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0) + PREVP=0. ! Rain evaporation (kg/kg; <0) + NSmICE=0. ! Cloud ice number concentration (m^-3) + Nrain=0. ! Rain number concentration (m^-3) !jul28 begin + Nsnow=0. ! "Snow" number concentration (m^-3) + RQRnew=0. ! Final rain content (kg/m**3) + RQSnew=0. ! Final "snow" content (kg/m**3) + Zrain=0. ! Radar reflectivity from rain (mm**6 m-3) + Zsnow=0. ! Radar reflectivity from snow (mm**6 m-3) + Ztot=0. ! Radar reflectivity from rain+snow (mm**6 m-3) + INDEXR=MDRmin ! Mean diameter of rain (microns) + INDEXR1=INDEXR ! 1st updated mean diameter of rain (microns) + INDEXR2=INDEXR ! 2nd updated mean diameter of rain (microns) + N0r=0. ! 1st estimate for rain intercept (m^-4) + DUM1=MIN(0.,TC) + DUM=XMImax*EXP(XMIexp*DUM1) + INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) ) ! 1st estimate for mean diameter of snow (microns) + VCI=0. ! Cloud ice fall speeds (m/s) + VSNOW=0. ! "Snow" (snow/graupel/sleet/hail) fall speeds (m/s) + VRAIN2=0. ! Rain fall speeds out of bottom of grid box (m/s) + RimeF1=1. ! Rime Factor (ratio, >=1, defined below) +! +!--- Double check input hydrometeor mixing ratios +! +! DUM=WC-(QI+QW+QR) +! DUM1=ABS(DUM) +! DUM2=TOLER*MIN(WC, QI+QW+QR) +! IF (DUM1 .GT. DUM2) THEN +! WRITE(0,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index, +! & ' L=',L +! WRITE(0,"(4(a12,g11.4,1x))") +! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC, +! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR +! ENDIF +! +!*********************************************************************** +!*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! **************** +!*********************************************************************** +! +!--- Calculate a few variables, which are used more than once below +! +!--- Latent heat of vaporization as a function of temperature from +! Bolton (1980, JAS) +! + TK2=1./(TK*TK) ! 1./TK**2 +! +!--- Basic thermodynamic quantities +! * DYNVIS - dynamic viscosity [ kg/(m*s) ] +! * THERM_COND - thermal conductivity [ J/(m*s*K) ] +! * DIFFUS - diffusivity of water vapor [ m**2/s ] +! + TFACTOR=SQRT(TK*TK*TK)/(TK+120.) + DYNVIS=1.496E-6*TFACTOR + THERM_COND=2.116E-3*TFACTOR + DIFFUS=8.794E-5*TK**1.81/PP +! +!--- Air resistance term for the fall speed of ice following the +! basic research by Heymsfield, Kajikawa, others +! + GAMMAS=MIN(1.5, (1.E5/PP)**C1) !-- limited to 1.5x +! +!--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470) +! + GAMMAR=(RHO0/RHO)**.4 +! +!---------------------------------------------------------------------- +!------------- IMPORTANT MICROPHYSICS DECISION TREE ----------------- +!---------------------------------------------------------------------- +! +!--- Determine if conditions supporting ice are present +! + IF (TC.LT.0. .OR. QI.GT. EPSQ .OR. ASNOW.GT.CLIMIT) THEN + ICE_logical=.TRUE. + ELSE + ICE_logical=.FALSE. + QLICE=0. + QTICE=0. + ENDIF + IF (T_ICE <= -100.) THEN + ICE_logical=.FALSE. + QLICE=0. + QTICE=0. + ENDIF +! +!--- Determine if rain is present +! + RAIN_logical=.FALSE. + IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE. +! +ice_test: IF (ICE_logical) THEN +! +!--- IMPORTANT: Estimate time-averaged properties. +! +!--- +! -> Small ice particles are assumed to have a mean diameter of 50 microns. +! * QSmICE - estimated mixing ratio for small cloud ice +!--- +! * TOT_ICE - total mass (small & large) ice before microphysics, +! which is the sum of the total mass of large ice in the +! current layer and the input flux of ice from above +! * PILOSS - greatest loss (<0) of total (small & large) ice by +! sublimation, removing all of the ice falling from above +! and the ice within the layer +! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed) +! ice mass to the unrimed ice mass (>=1) +! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1) +! * VSNOW - Fall speed of rimed snow w/ air resistance correction +! * VCI - Fall speed of 50-micron ice crystals w/ air resistance correction +! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer +! * XLIMASS - used for debugging, associated with calculating large ice mixing ratio +! * FLIMASS - mass fraction of large ice +! * QTICE - time-averaged mixing ratio of total ice +! * QLICE - time-averaged mixing ratio of large ice +! * NLICE - time-averaged number concentration of large ice +! * NSmICE - number concentration of small ice crystals at current level +! * QSmICE - mixing ratio of small ice crystals at current level +!--- +!--- Assumed number fraction of large ice particles to total (large & small) +! ice particles, which is based on a general impression of the literature. +! + NInuclei=0. + NSmICE=0. + QSmICE=0. + Rcw=0. + IF (TC<0.) THEN +! +!--- Max # conc of small ice crystals based on 10% of total ice content +! or the parameter NSI_max +! + NSImax=MAX(NSI_max, 0.1*RHO*QI/MASSI(MDImin) ) !aug27 +! +!-- Specify Fletcher, Cooper, Meyers, etc. here for ice nuclei concentrations +! Cooper (1986): NInuclei=MIN(5.*EXP(-0.304*TC), NSImax) +! Fletcher (1962): NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax) +! +!aug28: The formulas below mean that Fletcher is used for >-21C and Cooper at colder +! temperatures. In areas of high ice contents near the tops of deep convection, +! the number concentrations will be determined by the lower value of the "FQi" +! contribution to NSImax or Cooper. +! + NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax) !aug28 - Fletcher (1962) + NInuclei=MIN(5.*EXP(-0.304*TC), NInuclei) !aug28 - Cooper (1984) + IF (QI>EPSQ) THEN + DUM=RRHO*MASSI(MDImin) + NSmICE=MIN(NInuclei, QI/DUM) + QSmICE=NSmICE*DUM + ENDIF ! End IF (QI>EPSQ) + ENDIF ! End IF (TC<0.) + init_ice: IF (QI<=EPSQ .AND. ASNOW<=CLIMIT) THEN + TOT_ICE=0. + PILOSS=0. + RimeF1=1. + VrimeF=1. + VEL_INC=GAMMAS + VSNOW=0. + VSNOW1=0. + VCI=0. + EMAIRI=THICK + XLIMASS=RimeF1*MASSI(INDEXS) + FLIMASS=1. + QLICE=0. + RQLICE=0. + QTICE=0. + NLICE=0. + ELSE init_ice + ! + !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160), + ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships + ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66). + ! +! +!sep10 - Start of changes described in (23) at top of code. +! + TOT_ICE=THICK*QI+BLEND*ASNOW + PILOSS=-TOT_ICE/THICK + QLgICE=MAX(0., QI-QSmICE) !-- 1st-guess estimate of large ice + VCI=GAMMAS*VSNOWI(MDImin) +! +!-- Need to save this original value before two_pass iteration +! + LBEF=MAX(1,L-1) + RimeF1=(RimeF_col(L)*THICK*QLgICE & + & +RimeF_col(LBEF)*BLEND*ASNOW)/TOT_ICE +! +!mar08 see (32); !apr22a see (41) start - Estimate mean diameter (INDEXS in microns) + IF (RimeF1>2.) THEN + DUM3=RH_NgC*(RHO*QLgICE)**C1 !- convective mean diameter + DUM2=RH_NgT*(RHO*QLgICE)**C1 !- transition mean diameter + IF (RimeF1>=10.) THEN + DUM=DUM3 + ELSE IF (RimeF1>=5.) THEN + DUM=0.2*(RimeF1-5.) !- Blend at 5<=RF<10 + DUM=DUM3*DUM+DUM2*(1.-DUM) + ELSE + DUM1=REAL(INDEXS) !- stratiform mean diameter + DUM=0.33333*(RimeF1-2.) !- Blend at 2=5. .AND. INDEXS==MDImax .AND. RQLICE>RQhail) THEN +!- Additional increase using Thompson graupel/hail fall speeds + DUM=GAMMAS*AVhail*RQLICE**BVhail + IF (DUM>VSNOW) THEN + VSNOW=DUM + VEL_INC=VSNOW/VSNOWI(INDEXS) + ENDIF + ENDIF + XLIMASS=RimeF1*MASSI(INDEXS) + NLICE=RQLICE/XLIMASS +! +!sep16 - End of change described in (26) at top of code. +! +!=========================================== + IF (IPASS>=2 .OR. & + (NLICE>=NLImin .AND. NLICE<=NSI_max)) EXIT two_pass +!may17 - end +!=========================================== +! +!--- Force NLICE to be between NLImin and NSI_max when IPASS=1 +! +! IF (PRINT_diag .AND. RQLICE>Revised_LICE) DUM2=NLICE !-- For debugging (see DUM2 below) + NLICE=MAX(NLImin, MIN(NSI_max, NLICE) ) +!sep16 - End of changes +! + XLI=RQLICE/(NLICE*RimeF1) !- Mean mass of unrimed ice +new_size: IF (XLI<=MASSI(MDImin) ) THEN + INDEXS=MDImin + ELSE IF (XLI<=MASSI(450) ) THEN new_size + DLI=9.5885E5*XLI**.42066 ! DLI in microns + INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) ) + ELSE IF (XLIRevised_LICE) THEN +! WRITE(0,"(5(a12,g11.4,1x))") '{$ RimeF1=',RimeF1, & +! & ' RHO*QLICE=',RQLICE,' TC=',TC,' NLICE=',NLICE, & +! & ' NLICEold=',DUM2 +! Revised_LICE=1.2*RQLICE +! ENDIF + ENDIF new_size +!=========================================== + ENDDO two_pass +!=========================================== + ENDIF init_ice + ENDIF ice_test +! +!mar03 !may11 !jun01 - start; see (34) above + IF(TC<=0.) THEN + STRAT=.FALSE. + INDEXRmin=MDRmin + TIMLT=0. + INDEXS0C=INDEXS + RHOX0C=22.5*MAX(1.,MIN(RimeF1,40.)) !- Estimated ice density at 0C (kg m^-3) + ELSE ! TC>0. + IF(.NOT.RAIN_logical) THEN + STRAT=.FALSE. !- Reset STRAT + INDEXRmin=MDRmin !- Reset INDEXRmin + IF(.NOT.ICE_logical) TIMLT=0. + ELSE +!- STRAT=T for stratiform rain + IF(TIMLT>EPSQ .AND. RHOX0C<=225.) THEN + STRAT=.TRUE. + ELSE + STRAT=.FALSE. !- Reset STRAT + INDEXRmin=MDRmin !- Reset INDEXRmin + ENDIF + IF(STRAT .AND. INDEXRmin<=MDRmin) THEN + INDEXRmin=INDEXS0C*(0.001*RHOX0C)**C1 + INDEXRmin=MAX(MDRmin, MIN(INDEXRmin, INDEXRstrmax) ) + ENDIF + ENDIF + ENDIF +! + IF(STRAT .OR. TIMLT>EPSQ) THEN + DRZL=.FALSE. + ELSE +!- DRZL=T for drizzle (no melted ice falling from above) + DRZL=.TRUE. !mar30 + ENDIF +!jun01 - end +! +!---------------------------------------------------------------------- +!--------------- Calculate individual processes ----------------------- +!---------------------------------------------------------------------- +! +!--- Cloud water autoconversion to rain (PRAUT) and collection of cloud +! water by precipitation ice (PIACW) +! + IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN +!-- The old autoconversion threshold returns + DUM2=RHO*QW + IF (DUM2>QAUT0) THEN +!-- July 2010 version follows Liu & Daum (JAS, 2004) and Liu et al. (JAS, 2006) + DUM2=DUM2*DUM2 + DUM=BRAUT*DUM2*QW + DUM1=ARAUT*DUM2 + PRAUT=MIN(QW, DUM*(1.-EXP(-DUM1*DUM1)) ) + ENDIF + IF (QLICE .GT. EPSQ) THEN + ! + !--- Collection of cloud water by large ice particles ("snow") + ! PIACWI=PIACW for riming, PIACWI=0 for shedding + ! + FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS) ) !jul28 (16) + PIACW=FWS*QW + IF (TC<0.) THEN + PIACWI=PIACW !- Large ice riming + Rcw=ARcw*DUM2**C1 !- Cloud droplet radius in microns + ENDIF + ENDIF ! End IF (QLICE .GT. EPSQ) + ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) +! +!---------------------------------------------------------------------- +!--- Calculate homogeneous freezing of cloud water (PIACW, PIACWI) and +! ice deposition (PIDEP), which also includes ice initiation (PINIT) +! +ice_only: IF (TC.LT.T_ICE .AND. (WV.GT.QSWgrd .OR. QW.GT.EPSQ)) THEN + ! + !--- Adjust to ice saturation at T More extensive units conversion than can be described here to go from +! eq. (13) in Liu et al. (JAS, 2006) to what's programmed below. Note that +! the units used throughout the paper are in cgs units! +! + ARAUT=1.03e19/(NCW*SQRT(NCW)) + BRAUT=DTPH*1.1E10*BETA6/NCW +! +!--- QAUT0 is the *OLD* threshold cloud content for autoconversion to rain +! needed for droplets to reach a diameter of 20 microns (following +! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM). It's no longer +! used in this version, but the value is passed into radiation in case +! a ball park estimate is needed. +!--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations +! of 300, 200, and 100 cm**-3, respectively +! + QAUT0=PI*RHOL*NCW*(20.E-6)**3/6. !-- legacy +! +!--- For calculating cloud droplet radius in microns, Rcw +! + ARcw=1.E6*(0.75/(PI*NCW*RHOL))**C1 +! +!may20 - start +! +!- An explanation for the following settings assumed for "hail" or frozen drops (RF>10): +! RH_NgC=PI*500.*5.E3 +! RHOg=500 kg m^-3, Ng=5.e3 m^-3 => "hail" content >7.85 g m^-3 for INDEXS=MDImax +!- or - +! RH_NgC=PI*500.*1.E3 +! RHOg=500 kg m^-3, Ng=1.e3 m^-3 => "hail" content >1.57 g m^-3 for INDEXS=MDImax +! + RH_NgC=PI*500.*1.E3 !- RHOg=500 kg m^-3, Ng=1.e3 m^-3 + RQhail=RH_NgC*(1.E-3)**3 !- Threshold "hail" content ! becomes 1.57 g m^-3 + Bvhail=0.82*C1 !- Exponent for Thompson graupel + Avhail=1353.*(1./RH_NgC)**Bvhail !- 1353 ~ constant for Thompson graupel + RH_NgC=1.E6*(1./RH_NgC)**C1 !mar08 (convection, convert to microns) +! +!- An explanation for the following settings assumed for graupel (RF>5): +! RH_NgT=PI*300.*10.E3 +! RHOg=300 kg m^-3, Ng=10.e3 m^-3 => "graupel" content must exceed 9.43 g m^-3 for INDEXS=MDImax +!- or - +! RH_NgT=PI*300.*5.E3 +! RHOg=300 kg m^-3, Ng=5.e3 m^-3 => "graupel" content must exceed 4.71 g m^-3 for INDEXS=MDImax +! + RH_NgT=PI*300.*5.E3 !- RHOg=300 kg m^-3, Ng=5.e3 m^-3 + RH_NgT=1.E6*(1./RH_NgT)**C1 !mar08 (transition, convert to microns) +!may20 - end +! +!--- For calculating snow optical depths by considering bulk density of +! snow based on emails from Q. Fu (6/27-28/01), where optical +! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff +! is effective radius, and DENS is the bulk density of snow. +! +! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation +! T = 1.5*1.E3*SWPrad/(Reff*DENS) +! +! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW +! +! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3] +! + DO I=MDImin,MDImax + SDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I) + ENDDO +! + Thour_print=-DTPH/3600. +! + + RETURN +! +!----------------------------------------------------------------------- +! +9061 CONTINUE + WRITE(0,*)' module_mp_etanew: error opening ETAMPNEW_DATA.expanded_rain on unit ',etampnew_unit1 + STOP +! +!----------------------------------------------------------------------- + END SUBROUTINE FERRIER_INIT_hr +! +!>\ingroup hafs_famp + SUBROUTINE MY_GROWTH_RATES_NMM_hr (DTPH) +! +!--- Below are tabulated values for the predicted mass of ice crystals +! after 600 s of growth in water saturated conditions, based on +! calculations from Miller and Young (JAS, 1979). These values are +! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of +! Young (1993). Values at temperatures colder than -27C were +! assumed to be invariant with temperature. +! +!--- Used to normalize Miller & Young (1979) calculations of ice growth +! over large time steps using their tabulated values at 600 s. +! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993). +! + IMPLICIT NONE +! + REAL,INTENT(IN) :: DTPH +! + REAL DT_ICE + REAL,DIMENSION(35) :: MY_600 +!WRF +! +!----------------------------------------------------------------------- +!-- 20090714: These values are in g and need to be converted to kg below + DATA MY_600 / & + & 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, & + & 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, & + & 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, & + & 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, & + & 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, & + & 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, & + & 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C +! +!----------------------------------------------------------------------- +! + DT_ICE=(DTPH/600.)**1.5 + MY_GROWTH_NMM=DT_ICE*MY_600*1.E-3 !-- 20090714: Convert from g to kg +! +!----------------------------------------------------------------------- +! + END SUBROUTINE MY_GROWTH_RATES_NMM_hr +! +!----------------------------------------------------------------------- +!--------- Old GFS saturation vapor pressure lookup tables ----------- +!----------------------------------------------------------------------- +! +!>\ingroup hafs_famp + SUBROUTINE GPVS_hr +! ****************************************************************** +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: GPVS_hr COMPUTE SATURATION VAPOR PRESSURE TABLE +! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82 +! +! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF +! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS. +! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX. +! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH +! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN. +! +! PROGRAM HISTORY LOG: +! 91-05-07 IREDELL +! 94-12-30 IREDELL EXPAND TABLE +! 96-02-19 HONG ICE EFFECT +! 01-11-29 JIN MODIFIED FOR WRF +! +! USAGE: CALL GPVS_hr +! +! SUBPROGRAMS CALLED: +! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +! +!$$$ + IMPLICIT NONE + real :: X,XINC,T + integer :: JX +!---------------------------------------------------------------------- + XINC=(XMAX-XMIN)/(NX-1) + C1XPVS=1.-XMIN/XINC + C2XPVS=1./XINC + C1XPVS0=1.-XMIN/XINC + C2XPVS0=1./XINC +! + DO JX=1,NX + X=XMIN+(JX-1)*XINC + T=X + TBPVS(JX)=FPVSX(T) + TBPVS0(JX)=FPVSX0(T) + ENDDO +! + END SUBROUTINE GPVS_hr +!----------------------------------------------------------------------- +!*********************************************************************** +!----------------------------------------------------------------------- + REAL FUNCTION FPVS(T) +!----------------------------------------------------------------------- +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE +! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82 +! +! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE. +! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE +! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS. +! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA. +! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES. +! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION. +! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE. +! +! PROGRAM HISTORY LOG: +! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION +! 94-12-30 IREDELL EXPAND TABLE +! 96-02-19 HONG ICE EFFECT +! 01-11-29 JIN MODIFIED FOR WRF +! +! USAGE: PVS=FPVS(T) +! +! INPUT ARGUMENT LIST: +! T - REAL TEMPERATURE IN KELVIN +! +! OUTPUT ARGUMENT LIST: +! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB) +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +!$$$ + IMPLICIT NONE + real,INTENT(IN) :: T + real XJ + integer :: JX +!----------------------------------------------------------------------- + IF (T>=XMIN .AND. T<=XMAX) THEN + XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX)) + JX=MIN(XJ,NX-1.) + FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX)) + ELSE IF (T>XMAX) THEN +!-- Magnus Tetens formula for water saturation (Murray, 1967) +! (saturation vapor pressure in kPa) + FPVS=0.61078*exp(17.2694*(T-273.16)/(T-35.86)) + ELSE +!-- Magnus Tetens formula for ice saturation(Murray, 1967) +! (saturation vapor pressure in kPa) + FPVS=0.61078*exp(21.8746*(T-273.16)/(T-7.66)) + ENDIF +! + END FUNCTION FPVS +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + REAL FUNCTION FPVS0(T) +!----------------------------------------------------------------------- + IMPLICIT NONE + real,INTENT(IN) :: T + real :: XJ1 + integer :: JX1 +!----------------------------------------------------------------------- + IF (T>=XMIN .AND. T<=XMAX) THEN + XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX)) + JX1=MIN(XJ1,NX-1.) + FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1)) + ELSE +!-- Magnus Tetens formula for water saturation (Murray, 1967) +! (saturation vapor pressure in kPa) + FPVS0=0.61078*exp(17.2694*(T-273.16)/(T-35.86)) + ENDIF +! + END FUNCTION FPVS0 +!----------------------------------------------------------------------- +!*********************************************************************** +!----------------------------------------------------------------------- + REAL FUNCTION FPVSX(T) +!----------------------------------------------------------------------- +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE +! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82 +! +! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE. +! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS +! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID. +! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT +! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED. +! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT +! TO GET THE FORMULA +! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR)) +! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS +! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE. +! +! PROGRAM HISTORY LOG: +! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION +! 94-12-30 IREDELL EXACT COMPUTATION +! 96-02-19 HONG ICE EFFECT +! 01-11-29 JIN MODIFIED FOR WRF +! +! USAGE: PVS=FPVSX(T) +! REFERENCE: EMANUEL(1994),116-117 +! +! INPUT ARGUMENT LIST: +! T - REAL TEMPERATURE IN KELVIN +! +! OUTPUT ARGUMENT LIST: +! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB) +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +!$$$ + IMPLICIT NONE +!----------------------------------------------------------------------- + real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 & + , CLIQ=4.1855E+3,CVAP= 1.8460E+3 & + , CICE=2.1060E+3,HSUB=2.8340E+6 +! + real, parameter :: PSATK=PSAT*1.E-3 + real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP) + real, parameter :: DLDTI=CVAP-CICE & + , XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP) + real T,TR +!----------------------------------------------------------------------- + TR=TTP/T +! + IF(T.GE.TTP)THEN + FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR)) + ELSE + FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR)) + ENDIF +! + END FUNCTION FPVSX +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + REAL FUNCTION FPVSX0(T) +!----------------------------------------------------------------------- + IMPLICIT NONE + real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 & + , CLIQ=4.1855E+3,CVAP=1.8460E+3 & + , CICE=2.1060E+3,HSUB=2.8340E+6 + real,PARAMETER :: PSATK=PSAT*1.E-3 + real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP) + real,PARAMETER :: DLDTI=CVAP-CICE & + , XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP) + real :: T,TR +!----------------------------------------------------------------------- + TR=TTP/T + FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR)) +! + END FUNCTION FPVSX0 + +! + END MODULE module_mp_fer_hires diff --git a/physics/module_MYNNrad_post.F90 b/physics/module_MYNNrad_post.F90 index 7acd2e406..1364db62e 100644 --- a/physics/module_MYNNrad_post.F90 +++ b/physics/module_MYNNrad_post.F90 @@ -22,6 +22,7 @@ end subroutine mynnrad_post_finalize #endif SUBROUTINE mynnrad_post_run( & & ix,im,levs, & + & flag_init,flag_restart, & & qc,qi, & & qc_save, qi_save, & & errmsg, errflg ) @@ -34,6 +35,7 @@ SUBROUTINE mynnrad_post_run( & !------------------------------------------------------------------- integer, intent(in) :: ix, im, levs + logical, intent(in) :: flag_init, flag_restart real(kind=kind_phys), dimension(im,levs), intent(out) :: qc, qi real(kind=kind_phys), dimension(im,levs), intent(in) :: qc_save, qi_save character(len=*), intent(out) :: errmsg @@ -48,6 +50,11 @@ SUBROUTINE mynnrad_post_run( & !write(0,*)"==============================================" !write(0,*)"in mynn rad post" + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_post flag_init = ', flag_init + return + endif + ! Add subgrid cloud information: do k = 1, levs do i = 1, im diff --git a/physics/module_MYNNrad_post.meta b/physics/module_MYNNrad_post.meta index b09abe01e..881a19fff 100644 --- a/physics/module_MYNNrad_post.meta +++ b/physics/module_MYNNrad_post.meta @@ -25,6 +25,22 @@ type = integer intent = in optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F [qc] standard_name = cloud_condensed_water_mixing_ratio long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) diff --git a/physics/module_MYNNrad_pre.F90 b/physics/module_MYNNrad_pre.F90 index 858abebee..95dc95445 100644 --- a/physics/module_MYNNrad_pre.F90 +++ b/physics/module_MYNNrad_pre.F90 @@ -32,6 +32,7 @@ end subroutine mynnrad_pre_finalize !###=================================================================== SUBROUTINE mynnrad_pre_run( & & ix,im,levs, & + & flag_init,flag_restart, & & qc, qi, T3D, & & qc_save, qi_save, & & qc_bl,cldfra_bl, & @@ -50,6 +51,7 @@ SUBROUTINE mynnrad_pre_run( & ! Interface variables real (kind=kind_phys), parameter :: gfac=1.0e5/con_g integer, intent(in) :: ix, im, levs + logical, intent(in) :: flag_init, flag_restart real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc, qi real(kind=kind_phys), dimension(im,levs), intent(in) :: T3D,delp real(kind=kind_phys), dimension(im,levs), intent(inout) :: & @@ -71,13 +73,17 @@ SUBROUTINE mynnrad_pre_run( & !write(0,*)"==============================================" !write(0,*)"in mynn rad pre" + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_pre flag_init = ', flag_init + return + endif ! Add subgrid cloud information: do k = 1, levs do i = 1, im qc_save(i,k) = qc(i,k) qi_save(i,k) = qi(i,k) - clouds1(i,k)=CLDFRA_BL(i,k) + clouds1(i,k) = CLDFRA_BL(i,k) IF (qc(i,k) < 1.E-6 .AND. qi(i,k) < 1.E-8 .AND. CLDFRA_BL(i,k)>0.001) THEN !Partition the BL clouds into water & ice according to a linear diff --git a/physics/module_MYNNrad_pre.meta b/physics/module_MYNNrad_pre.meta index 617ee3f31..3b5943c66 100644 --- a/physics/module_MYNNrad_pre.meta +++ b/physics/module_MYNNrad_pre.meta @@ -25,6 +25,22 @@ type = integer intent = in optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F [qc] standard_name = cloud_condensed_water_mixing_ratio long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index ea5800736..7345f2667 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -67,7 +67,7 @@ SUBROUTINE LSMRUC( & Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & GLW,GSW,EMISS,CHKLOWQ, CHS, & FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, & - Z0,SNOALB,ALBBCK, & !Z0,SNOALB,ALBBCK,LAI, & + Z0,SNOALB,ALBBCK,LAI, & landusef, nlcat, & ! mosaic_lu, mosaic_soil, & soilctop, nscat, & QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, & @@ -218,6 +218,7 @@ SUBROUTINE LSMRUC( & CANWAT, & ! new SNOALB, & ALB, & + LAI, & EMISS, & MAVAIL, & SFCEXC, & @@ -269,7 +270,6 @@ SUBROUTINE LSMRUC( & PC, & SFCRUNOFF, & UDRUNOFF, & - LAI, & EMISSL, & ZNTL, & LMAVAIL, & @@ -431,8 +431,8 @@ SUBROUTINE LSMRUC( & !! or ~100 mm of snow height ! ! snowc(i,j) = min(1.,snow(i,j)/32.) - soilt1(i,j)=soilt(i,j) - if(snow(i,j).le.32.) soilt1(i,j)=tso(i,1,j) +! soilt1(i,j)=soilt(i,j) +! if(snow(i,j).le.32.) soilt1(i,j)=tso(i,1,j) !> - Initializing inside snow temp if it is not defined IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN IF(snow(i,j).gt.32.) THEN @@ -450,7 +450,9 @@ SUBROUTINE LSMRUC( & patmb=P8w(i,kms,j)*1.e-2 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATMB IF((qvg(i,j) .LE. 0.) .or. (qvg(i,j) .GT.0.1)) THEN - qvg (i,j) = QSG(i,j)*mavail(i,j) + !17sept19 - bad approximation with very low mavail. + !qvg(i,j) = QSG(i,j)*mavail(i,j) + qvg (i,j) = qv3d(i,1,j) IF (debug_print ) THEN print *, & 'QVG is initialized in RUCLSM ', qvg(i,j),mavail(i,j),qsg(i,j),i,j @@ -751,7 +753,7 @@ SUBROUTINE LSMRUC( & meltfactor = 0.85 do k=2,nzs - if(zsmain(k).ge.1.0) then + if(zsmain(k).ge.1.1) then NROOT=K goto 111 endif diff --git a/physics/mp_fer_hires.F90 b/physics/mp_fer_hires.F90 new file mode 100644 index 000000000..9f265db22 --- /dev/null +++ b/physics/mp_fer_hires.F90 @@ -0,0 +1,401 @@ +!>\file mp_fer_hires.F90 +!! This file contains + +! +module mp_fer_hires + + use machine, only : kind_phys + + use module_mp_fer_hires, only : ferrier_init_hr, FER_HIRES + + implicit none + + public :: mp_fer_hires_init, mp_fer_hires_run, mp_fer_hires_finalize + + private + + logical :: is_initialized = .False. + + ! * T_ICE - temperature (C) threshold at which all remaining liquid water + ! is glaciated to ice + ! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs + REAL, PUBLIC, PARAMETER :: T_ICE=-40., & + T0C=273.15, & + T_ICEK=T0C+T_ICE + + contains + +!> This subroutine initialize constants & lookup tables for Ferrier-Aligao MP +!! scheme. +!> \section arg_table_mp_fer_hires_init Argument Table +!! \htmlinclude mp_fer_hires_init.html +!! + subroutine mp_fer_hires_init(ncol, nlev, dtp, imp_physics, & + imp_physics_fer_hires, & + restart, & + f_ice,f_rain,f_rimef, & + mpicomm, mpirank,mpiroot, & + threads, errmsg, errflg) + + USE machine, ONLY : kind_phys + USE MODULE_MP_FER_HIRES, ONLY : FERRIER_INIT_HR + implicit none + + integer, intent(in) :: ncol + integer, intent(in) :: nlev + real(kind_phys), intent(in) :: dtp + integer, intent(in) :: imp_physics + integer, intent(in) :: imp_physics_fer_hires + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + integer, intent(in) :: threads + logical, intent(in) :: restart + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + real(kind_phys), intent(out), optional :: f_ice(1:ncol,1:nlev) + real(kind_phys), intent(out), optional :: f_rain(1:ncol,1:nlev) + real(kind_phys), intent(out), optional :: f_rimef(1:ncol,1:nlev) + + + ! Local variables + integer :: ims, ime, lm,i,k + !real(kind=kind_phys) :: DT_MICRO + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + ! Set internal dimensions + ims = 1 + ime = ncol + lm = nlev + + ! MZ* temporary + if (mpirank==mpiroot) then + write(0,*) ' -----------------------------------------------' + write(0,*) ' --- !!! WARNING !!! ---' + write(0,*) ' --- the CCPP Ferrier-Aligo MP scheme is ---' + write(0,*) ' --- currently under development, use at ---' + write(0,*) ' --- your own risk . ---' + write(0,*) ' -----------------------------------------------' + end if + ! MZ* temporary + + if (imp_physics /= imp_physics_fer_hires ) then + write(errmsg,'(*(a))') "Logic error: namelist choice of microphysics is different from Ferrier-Aligo MP" + errflg = 1 + return + end if + + !MZ: fer_hires_init() in HWRF + IF(.NOT.RESTART .AND. present(F_ICE)) THEN !HWRF + write(errmsg,'(*(a))') " WARNING: F_ICE,F_RAIN AND F_RIMEF IS REINITIALIZED " + DO K = 1,lm + DO I= ims,ime + F_ICE(i,k)=0. + F_RAIN(i,k)=0. + F_RIMEF(i,k)=1. + ENDDO + ENDDO + ENDIF + !MZ: fer_hires_init() in HWRF + + CALL FERRIER_INIT_HR(dtp,mpicomm,mpirank,mpiroot,threads) + + if (mpirank==mpiroot) write (0,*)'F-A: FERRIER_INIT_HR finished ...' + if (errflg /= 0 ) return + + is_initialized = .true. + + + end subroutine mp_fer_hires_init + +!>\defgroup hafs_famp HAFS Ferrier-Aligo Cloud Microphysics Scheme +!> This is the CCPP-compliant FER_HIRES driver module. +!> \section arg_table_mp_fer_hires_run Argument Table +!! \htmlinclude mp_fer_hires_run.html +!! + SUBROUTINE mp_fer_hires_run(NCOL, NLEV, DT ,SPEC_ADV & + ,SLMSK & + ,PRSI,P_PHY & + ,T,Q,CWM & + ,TRAIN,SR & + ,F_ICE,F_RAIN,F_RIMEF & + ,QC,QR,QI,QG & ! wet mixing ratio + !,qc_m,qi_m,qr_m & + ,PREC &!,ACPREC -MZ:not used + ,mpirank, mpiroot, threads & + ,refl_10cm & + ,RHGRD,dx & + ,EPSQ,R_D,P608,CP,G & + ,errmsg,errflg) + +!----------------------------------------------------------------------- + USE MACHINE, ONLY: kind_phys +! + IMPLICIT NONE +! +!----------------------------------------------------------------------- +! + INTEGER,PARAMETER :: D_SS=1 +! +!------------------------ +!*** Argument Variables +!------------------------ + + integer, intent(in ) :: ncol + integer, intent(in ) :: nlev + real(kind_phys), intent(in ) :: dt + integer, intent(in ) :: threads + logical, intent(in ) :: spec_adv + integer, intent(in ) :: mpirank + integer, intent(in ) :: mpiroot + real(kind_phys), intent(in ) :: slmsk(1:ncol) + real(kind_phys), intent(in ) :: prsi(1:ncol,1:nlev+1) + real(kind_phys), intent(in ) :: p_phy(1:ncol,1:nlev) + real(kind_phys), intent(in ) :: epsq,r_d,p608,cp,g + real(kind_phys), intent(inout) :: t(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: q(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: cwm(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: train(1:ncol,1:nlev) + real(kind_phys), intent(out ) :: sr(1:ncol) + real(kind_phys), intent(inout) :: f_ice(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: f_rain(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: f_rimef(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qc(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qr(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qi(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qg(1:ncol,1:nlev) ! QRIMEF + + real(kind_phys), intent(inout) :: prec(1:ncol) +! real(kind_phys) :: acprec(1:ncol) !MZ: change to local + real(kind_phys), intent(inout) :: refl_10cm(1:ncol,1:nlev) + real(kind_phys), intent(in ) :: rhgrd + real(kind_phys), intent(in ) :: dx(1:ncol) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg +! +!--------------------- +!*** Local Variables +!--------------------- +! + integer :: I,J,K,N + integer :: lowlyr(1:ncol) + integer :: dx1 + !real(kind_phys) :: mprates(1:ncol,1:nlev,d_ss) + real(kind_phys) :: DTPHS,PCPCOL,RDTPHS,TNEW + real(kind_phys) :: ql(1:nlev),tl(1:nlev) + real(kind_phys) :: rainnc(1:ncol),rainncv(1:ncol) + real(kind_phys) :: snownc(1:ncol),snowncv(1:ncol) + real(kind_phys) :: graupelncv(1:ncol) + real(kind_phys) :: dz(1:ncol,1:nlev) + real(kind_phys) :: pi_phy(1:ncol,1:nlev) + real(kind_phys) :: rr(1:ncol,1:nlev) + real(kind_phys) :: th_phy(1:ncol,1:nlev) + real(kind_phys) :: R_G, CAPPA + +! Dimension + integer :: ims, ime, jms, jme, lm + +!----------------------------------------------------------------------- +!*********************************************************************** +!----------------------------------------------------------------------- + R_G=1./G + CAPPA=R_D/CP + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Check initialization state + if (.not. is_initialized) then + write(errmsg, fmt='((a))') 'mp_fer_hires_run called before mp_fer_hires_init' + errflg = 1 + return + end if + + +!ZM NTSD=ITIMESTEP +!ZM presume nphs=1 DTPHS=NPHS*DT + DTPHS=DT + RDTPHS=1./DTPHS +!ZM AVRAIN=AVRAIN+1. + +! Set internal dimensions + ims = 1 + ime = ncol + jms = 1 + jme = 1 + lm = nlev + +! Use the dx of the 1st i point to set an integer value of dx to be used for +! determining where RHgrd should be set to 0.98 in the coarse domain when running HAFS. + DX1=NINT(DX(1)) + +!----------------------------------------------------------------------- +!*** NOTE: THE NMMB HAS IJK STORAGE WITH LAYER 1 AT THE TOP. +!*** THE WRF PHYSICS DRIVERS HAVE IKJ STORAGE WITH LAYER 1 +!*** AT THE BOTTOM. +!----------------------------------------------------------------------- +!....................................................................... + DO I=IMS,IME +! + LOWLYR(I)=1 +! +!----------------------------------------------------------------------- +!*** FILL RAINNC WITH ZERO (NORMALLY CONTAINS THE NONCONVECTIVE +!*** ACCUMULATED RAIN BUT NOT YET USED BY NMM) +!*** COULD BE OBTAINED FROM ACPREC AND CUPREC (ACPREC-CUPREC) +!----------------------------------------------------------------------- +!..The NC variables were designed to hold simulation total accumulations +!.. whereas the NCV variables hold timestep only values, so change below +!.. to zero out only the timestep amount preparing to go into each +!.. micro routine while allowing NC vars to accumulate continually. +!.. But, the fact is, the total accum variables are local, never saved +!.. nor written so they go nowhere at the moment. +! + RAINNC (I)=0. ! NOT YET USED BY NMM + RAINNCv(I)=0. + SNOWNCv(I)=0. + graupelncv(i) = 0.0 +! +!----------------------------------------------------------------------- +!*** FILL THE SINGLE-COLUMN INPUT +!----------------------------------------------------------------------- +! + DO K=LM,1,-1 ! We are moving down from the top in the flipped arrays + +! +! TL(K)=T(I,K) +! QL(K)=AMAX1(Q(I,K),EPSQ) +! + RR(I,K)=P_PHY(I,K)/(R_D*T(I,K)*(P608*AMAX1(Q(I,K),EPSQ)+1.)) + PI_PHY(I,K)=(P_PHY(I,K)*1.E-5)**CAPPA + TH_PHY(I,K)=T(I,K)/PI_PHY(I,K) + DZ(I,K)=(PRSI(I,K)-PRSI(I,K+1))*R_G/RR(I,K) + +! +!*** CALL MICROPHYSICS + +!MZ* in HWRF +!-- 6/11/2010: Update cwm, F_ice, F_rain and F_rimef arrays + cwm(I,K)=QC(I,K)+QR(I,K)+QI(I,K) + IF (QI(I,K) <= EPSQ) THEN + F_ICE(I,K)=0. + F_RIMEF(I,K)=1. + IF (T(I,K) < T_ICEK) F_ICE(I,K)=1. + ELSE + F_ICE(I,K)=MAX( 0., MIN(1., QI(I,K)/cwm(I,K) ) ) + F_RIMEF(I,K)=QG(I,K)/QI(I,K) + ENDIF + IF (QR(I,K) <= EPSQ) THEN + F_RAIN(I,K)=0. + ELSE + F_RAIN(I,K)=QR(I,K)/(QR(I,K)+QC(I,K)) + ENDIF + + end do + enddo + +!--------------------------------------------------------------------- +!*** Update the rime factor array after 3d advection +!--------------------------------------------------------------------- +!MZ* in namphysics +! DO K=1,LM +! DO I=IMS,IME +! IF (QG(I,K)>EPSQ .AND. QI(I,K)>EPSQ) THEN +! F_RIMEF(I,K)=MIN(50.,MAX(1.,QG(I,K)/QI(I,K))) +! ELSE +! F_RIMEF(I,K)=1. +! ENDIF +! ENDDO +! ENDDO + + +!--------------------------------------------------------------------- + + CALL FER_HIRES( & + DT=dtphs,RHgrd=RHGRD & + ,DZ8W=dz,RHO_PHY=rr,P_PHY=p_phy,PI_PHY=pi_phy & + ,TH_PHY=th_phy,T_PHY=t & + ,Q=Q,QT=cwm & + ,LOWLYR=LOWLYR,SR=SR & + ,F_ICE_PHY=F_ICE,F_RAIN_PHY=F_RAIN & + ,F_RIMEF_PHY=F_RIMEF & + ,QC=QC,QR=QR,QS=QI & + ,RAINNC=rainnc,RAINNCV=rainncv & + ,threads=threads & + ,IMS=IMS,IME=IME,JMS=JMS,JME=JME,LM=LM & + ,D_SS=d_ss & + ,refl_10cm=refl_10cm,DX1=DX1) + + +!....................................................................... + +!MZ* +!Aligo Oct-23-2019 +! - Convert dry qc,qr,qi back to wet mixing ratio +! DO K = 1, LM +! DO I= IMS, IME +! qc_m(i,k) = qc(i,k)/(1.0_kind_phys+q(i,k)) +! qi_m(i,k) = qi(i,k)/(1.0_kind_phys+q(i,k)) +! qr_m(i,k) = qr(i,k)/(1.0_kind_phys+q(i,k)) +! ENDDO +! ENDDO + + + +!----------------------------------------------------------- + DO K=1,LM + DO I=IMS,IME + +!--------------------------------------------------------------------- +!*** Calculate graupel from total ice array and rime factor +!--------------------------------------------------------------------- + +!MZ + IF (SPEC_ADV) then + QG(I,K)=QI(I,K)*F_RIMEF(I,K) + ENDIF + +! +!----------------------------------------------------------------------- +!*** UPDATE TEMPERATURE, SPECIFIC HUMIDITY, CLOUD WATER, AND HEATING. +!----------------------------------------------------------------------- +! + TNEW=TH_PHY(I,K)*PI_PHY(I,K) + TRAIN(I,K)=TRAIN(I,K)+(TNEW-T(I,K))*RDTPHS + T(I,K)=TNEW + ENDDO + ENDDO + +!....................................................................... + +! +!----------------------------------------------------------------------- +!*** UPDATE PRECIPITATION +!----------------------------------------------------------------------- +! + DO I=IMS,IME + PCPCOL=RAINNCV(I)*1.E-3 !MZ:unit:m + PREC(I)=PREC(I)+PCPCOL +!MZ ACPREC(I)=ACPREC(I)+PCPCOL !MZ: not used +! +! NOTE: RAINNC IS ACCUMULATED INSIDE MICROPHYSICS BUT NMM ZEROES IT OUT ABOVE +! SINCE IT IS ONLY A LOCAL ARRAY FOR NOW +! + ENDDO +!----------------------------------------------------------------------- +! + end subroutine mp_fer_hires_run + + +!> \section arg_table_mp_fer_hires_finalize Argument Table +!! + subroutine mp_fer_hires_finalize () + end subroutine mp_fer_hires_finalize + +end module mp_fer_hires diff --git a/physics/mp_fer_hires.meta b/physics/mp_fer_hires.meta new file mode 100644 index 000000000..36b40a95c --- /dev/null +++ b/physics/mp_fer_hires.meta @@ -0,0 +1,426 @@ +[ccpp-arg-table] + name = mp_fer_hires_init + type = scheme +[ncol] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[nlev] + standard_name = vertical_dimension + long_name = vertical layer dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[dtp] + standard_name = time_step_for_physics + long_name = physics timestep + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[imp_physics] + standard_name = flag_for_microphysics_scheme + long_name = choice of microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imp_physics_fer_hires] + standard_name = flag_for_fer_hires_microphysics_scheme + long_name = choice of Ferrier-Aligo microphysics scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F +[f_ice] + standard_name = fraction_of_ice_water_cloud + long_name = fraction of ice water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = T +[f_rain] + standard_name = fraction_of_rain_water_cloud + long_name = fraction of rain water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = T +[f_rimef] + standard_name = rime_factor + long_name = rime factor + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = T +[mpicomm] + standard_name = mpi_comm + long_name = MPI communicator + units = index + dimensions = () + type = integer + intent = in + optional = F +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in + optional = F +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in + optional = F +[threads] + standard_name = omp_threads + long_name = number of OpenMP threads available to scheme + units = count + dimensions = () + type = integer + intent = in + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F +######################################################################## +[ccpp-arg-table] + name = mp_fer_hires_finalize + type = scheme +######################################################################## +[ccpp-arg-table] + name = mp_fer_hires_run + type = scheme +[ncol] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[nlev] + standard_name = vertical_dimension + long_name = vertical layer dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[dt] + standard_name = time_step_for_physics + long_name = physics time step + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[spec_adv] + standard_name = flag_for_individual_cloud_species_advected + long_name = flag for individual cloud species advected + units = flag + dimensions = () + type = logical + intent = in + optional = F +[slmsk] + standard_name = sea_land_ice_mask_real + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_dimension) + type = real + kind= kind_phys + intent = in + optional = F +[prsi] + standard_name = air_pressure_at_interface + long_name = air pressure at model layer interfaces + units = Pa + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[p_phy] + standard_name = air_pressure + long_name = mean layer pressure + units = Pa + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[t] + standard_name = air_temperature_updated_by_physics + long_name = temperature updated by physics + units = K + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[q] + standard_name = water_vapor_specific_humidity_updated_by_physics + long_name = water vapor specific humidity updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cwm] + standard_name = total_cloud_condensate_mixing_ratio_updated_by_physics + long_name = total cloud condensate mixing ratio (except water vapor) updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[train] + standard_name = accumulated_change_of_air_temperature_due_to_FA_scheme + long_name = accumulated change of air temperature due to FA MP scheme + units = K + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sr] + standard_name = ratio_of_snowfall_to_rainfall + long_name = snow ratio: ratio of snow to total precipitation (explicit only) + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[f_ice] + standard_name = fraction_of_ice_water_cloud + long_name = fraction of ice water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[f_rain] + standard_name = fraction_of_rain_water_cloud + long_name = fraction of rain water cloud + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[f_rimef] + standard_name = rime_factor + long_name = rime factor + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[qc] + standard_name = cloud_condensed_water_mixing_ratio_updated_by_physics + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud condensed water updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qi] + standard_name = ice_water_mixing_ratio_updated_by_physics + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qr] + standard_name = rain_water_mixing_ratio_updated_by_physics + long_name = moist (dry+vapor, no condensates) mixing ratio of rain water updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qg] + standard_name = mass_weighted_rime_factor_updated_by_physics + long_name = mass weighted rime factor updated by physics + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[prec] + standard_name = lwe_thickness_of_explicit_precipitation_amount + long_name = explicit precipitation ( rain, ice, snow, graupel, ...) on physics timestep + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in + optional = F +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in + optional = F +[threads] + standard_name = omp_threads + long_name = number of OpenMP threads available to scheme + units = count + dimensions = () + type = integer + intent = in + optional = F +[refl_10cm] + standard_name = radar_reflectivity_10cm + long_name = instantaneous refl_10cm + units = dBZ + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rhgrd] + standard_name = fa_threshold_relative_humidity_for_onset_of_condensation + long_name = relative humidity threshold parameter for condensation for FA scheme + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[dx] + standard_name = cell_size + long_name = relative dx for the grid cell + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[EPSQ] + standard_name = minimum_value_of_specific_humidity + long_name = floor value for specific humidity + units = kg kg-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[R_D] + standard_name = gas_constant_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[P608] + standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one + long_name = (rv/rd) - 1 (rv = ideal gas constant for water vapor) + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[CP] + standard_name = specific_heat_of_dry_air_at_constant_pressure + long_name = specific heat of dry air at constant pressure + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[G] + standard_name = gravitational_acceleration + long_name = gravitational acceleration + units = m s-2 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 258c3412e..189b28c4c 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -262,6 +262,7 @@ module module_radiation_clouds !!\n =8: Thompson microphysics !!\n =6: WSM6 microphysics !!\n =10: MG microphysics +!!\n =15: Ferrier-Aligo microphysics !!\param me print control flag !>\section gen_cld_init cld_init General Algorithm !! @{ @@ -350,6 +351,8 @@ subroutine cld_init & print *,' --- WSM6 cloud microphysics' elseif (imp_physics == 10) then print *,' --- MG cloud microphysics' + elseif (imp_physics == 15) then + print *,' --- Ferrier-Aligo cloud microphysics' else print *,' !!! ERROR in cloud microphysc specification!!!', & & ' imp_physics (NP3D) =',imp_physics diff --git a/physics/sfc_drv_ruc.F90 b/physics/sfc_drv_ruc.F90 index fe12b5e17..3b4b8a118 100644 --- a/physics/sfc_drv_ruc.F90 +++ b/physics/sfc_drv_ruc.F90 @@ -139,11 +139,11 @@ end subroutine lsm_ruc_finalize ! DH* TODO - make order of arguments the same as in the metadata table subroutine lsm_ruc_run & ! inputs & ( iter, me, master, kdt, im, nlev, lsoil_ruc, lsoil, zs, & - & t1, q1, qc, soiltyp, vegtype, sigmaf, & + & t1, q1, qc, soiltyp, vegtype, sigmaf, laixy, & & sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, & & prsl1, zf, wind, shdmin, shdmax, alvwf, alnwf, & & snoalb, sfalb, flag_iter, flag_guess, isot, ivegsrc, fice, & - & smc, stc, slc, lsm_ruc, lsm, land, islimsk, & + & smc, stc, slc, lsm_ruc, lsm, land, islimsk, rdlai, & & imp_physics, imp_physics_gfdl, imp_physics_thompson, & & smcwlt2, smcref2, do_mynnsfclay, & & con_cp, con_rv, con_rd, con_g, con_pi, con_hvap, con_fvirt,& ! constants @@ -178,6 +178,8 @@ subroutine lsm_ruc_run & ! inputs & ch, prsl1, wind, shdmin, shdmax, & & snoalb, alvwf, alnwf, zf, qc, q1 + real (kind=kind_phys), dimension(:), intent(in) :: laixy + real (kind=kind_phys), intent(in) :: delt real (kind=kind_phys), intent(in) :: con_cp, con_rv, con_g, & con_pi, con_rd, & @@ -187,6 +189,8 @@ subroutine lsm_ruc_run & ! inputs integer, dimension(im), intent(in) :: islimsk ! sea/land/ice mask (=0/1/2) logical, intent(in) :: do_mynnsfclay + logical, intent(in) :: rdlai + ! --- in/out: integer, dimension(im), intent(inout) :: soiltyp, vegtype real (kind=kind_phys), dimension(lsoil_ruc) :: dzs @@ -317,6 +321,8 @@ subroutine lsm_ruc_run & ! inputs zs, sh2o, smfrkeep, tslb, smois, wetness, & ! out me, master, errmsg, errflg) + xlai = 0. + endif ! flag_init=.true.,iter=1 !-- end of initialization @@ -508,10 +514,10 @@ subroutine lsm_ruc_run & ! inputs ffrozp(i,j) = real(nint(srflag(i)),kind_phys) endif - !tgs - for now set rdlai2d to .false., WRF has LAI maps, and RUC LSM - ! uses rdlai2d = .true. - rdlai2d = .false. - !if( .not. rdlai2d) xlai = lai_data(vtype) + !tgs - rdlai is .false. when the LAI data is not available in the + ! - INPUT/sfc_data.nc + + rdlai2d = rdlai conflx2(i,1,j) = zf(i) * 2. ! factor 2. is needed to get the height of ! atm. forcing inside RUC LSM (inherited @@ -552,13 +558,15 @@ subroutine lsm_ruc_run & ! inputs !prcp(i,j) = rhoh2o * tprcp(i) ! tprcp in [m] - convective plus explicit !raincv(i,j) = rhoh2o * rainc(i) ! total time-step convective precip !rainncv(i,j) = rhoh2o * max(rain(i)-rainc(i),0.0) ! total time-step explicit precip - !graupelncv(i,j) = rhoh2o * graupel(i) - !snowncv(i,j) = rhoh2o * snow(i) - prcp(i,j) = rhoh2o * (rainc(i)+rainnc(i)) ! tprcp in [m] - convective plus explicit - raincv(i,j) = rhoh2o * rainc(i) ! total time-step convective precip - rainncv(i,j) = rhoh2o * rainnc(i) ! total time-step explicit precip + prcp(i,j) = rhoh2o * (rainc(i)+rainnc(i)) ! [mm] - convective plus explicit + raincv(i,j) = rhoh2o * rainc(i) ! [mm] - total time-step convective precip + rainncv(i,j) = rhoh2o * rainnc(i) ! [mm] - total time-step explicit precip graupelncv(i,j) = rhoh2o * graupel(i) snowncv(i,j) = rhoh2o * snow(i) + !if(prcp(i,j) > 0. .and. i==21) then + !print *,'prcp(i,j),rainncv(i,j),graupelncv(i,j),snowncv(i,j),ffrozp(i,j)',i,j, & + ! prcp(i,j),rainncv(i,j),graupelncv(i,j),snowncv(i,j),ffrozp(i,j) + !endif ! ice not used ! precipfr(i,j) = rainncv(i,j) * ffrozp(i,j) @@ -601,6 +609,8 @@ subroutine lsm_ruc_run & ! inputs albbck(i,j) = max(0.01, 0.5 * (alvwf(i) + alnwf(i))) alb(i,j) = sfalb(i) + if(rdlai2d) xlai(i,j) = laixy(i) + tbot(i,j) = tg3(i) !> - 4. history (state) variables (h): @@ -686,7 +696,7 @@ subroutine lsm_ruc_run & ! inputs znt(i,j) = zorl(i)/100. if(debug_print) then - !if(me==0 .and. i==ipr) then + if(me==0 .and. i==ipr) then write (0,*)'before RUC smsoil = ',smsoil(i,:,j), i,j write (0,*)'stsoil = ',stsoil(i,:,j), i,j write (0,*)'soilt = ',soilt(i,j), i,j @@ -780,7 +790,7 @@ subroutine lsm_ruc_run & ! inputs write (0,*)'shdmin1d(i,j) =',i,j,shdmin1d(i,j) write (0,*)'shdmax1d(i,j) =',i,j,shdmax1d(i,j) write (0,*)'rdlai2d =',rdlai2d - !endif + endif endif !> - Call RUC LSM lsmruc(). @@ -796,8 +806,7 @@ subroutine lsm_ruc_run & ! inputs & chs(i,j), flqc(i,j), flhc(i,j), & ! --- input/outputs: & wet(i,j), cmc(i,j), shdfac(i,j), alb(i,j), znt(i,j), & - & z0(i,j), snoalb1d(i,j), albbck(i,j), & -! & z0, snoalb1d, alb, xlai, & + & z0(i,j), snoalb1d(i,j), albbck(i,j), xlai(i,j), & & landusef(i,:,j), nlcat, & ! --- mosaic_lu and mosaic_soil are moved to the namelist ! & mosaic_lu, mosaic_soil, & diff --git a/physics/sfc_drv_ruc.meta b/physics/sfc_drv_ruc.meta index dac459405..3ae9a57a3 100644 --- a/physics/sfc_drv_ruc.meta +++ b/physics/sfc_drv_ruc.meta @@ -198,6 +198,12 @@ type = integer intent = in optional = F +[rdlai] + standard_name = flag_for_reading_leaf_area_index_from_input + long_name = flag for reading leaf area index from initial conditions for RUC LSM + units = flag + dimensions = () + type = logical [zs] standard_name = depth_of_soil_levels_for_land_surface_model long_name = depth of soil levels for land surface model @@ -529,6 +535,14 @@ kind = kind_phys intent = in optional = F +[laixy] + standard_name = leaf_area_index + long_name = leaf area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + optional = F [sfalb] standard_name = surface_diffused_shortwave_albedo long_name = mean surface diffused sw albedo diff --git a/physics/ugwp_driver_v0.F b/physics/ugwp_driver_v0.F index 52375dd18..7f5490d24 100644 --- a/physics/ugwp_driver_v0.F +++ b/physics/ugwp_driver_v0.F @@ -1993,8 +1993,8 @@ subroutine edmix_ugwp_v0(im, levs, dtp, Km(1:levs) = ksum(1:levs) * rho(1:levs)* rho(1:levs) do j=1, nstab - call diff_1d_wtend(levs, dtstab, Fw, Fw1, levs, - & del(i,:), Sw, Sw1) + call diff_1d_wtend(levs, dtstab, Fw, Fw1, Km, + & rdp, rdpm, Sw, Sw1) Fw = Sw Fw1 = Sw1 enddo @@ -2006,7 +2006,7 @@ subroutine edmix_ugwp_v0(im, levs, dtp, Kpt = Km*iPr_pt Fw(1:levs) = pdTdt(i, 1:levs)*Ptmap(1:levs) do j=1, nstab - call diff_1d_ptend(levs, dtstab, Fw, Kpt, del(i,:), Sw) + call diff_1d_ptend(levs, dtstab, Fw, Kpt, rdp, rdpm, Sw) Fw = Sw enddo ed_dtdt(i,1:levs) = Sw(1:levs)/Ptmap(1:levs)