diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index b9f0ea8bc..b4fe70027 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -2537,10 +2537,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c c------- final changed variable per unit mass flux c -!> - If grid size is less than a threshold value (dxcrtas: currently 8km if progsigma is not used and 30km if progsigma is used), the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. +!> - If grid size is less than a threshold value (dxcrtas: currently 8km if progsigma is not used), or progsigma = true, the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. ! if(progsigma)then - dxcrtas=30.e3 + dxcrtas=500.e3 dxcrtuf=10.e3 else dxcrtas=8.e3 diff --git a/physics/CONV/SAMF/samfdeepcnv.meta b/physics/CONV/SAMF/samfdeepcnv.meta index 7c208d2d3..e6f44bc7b 100644 --- a/physics/CONV/SAMF/samfdeepcnv.meta +++ b/physics/CONV/SAMF/samfdeepcnv.meta @@ -472,7 +472,7 @@ [omegain] standard_name = prognostic_updraft_velocity_in_convection long_name = convective updraft velocity - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys @@ -481,7 +481,7 @@ [omegaout] standard_name = updraft_velocity_updated_by_physics long_name = convective updraft velocity updated by physics - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 0746d085e..ee7f6a76e 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -204,7 +204,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & parameter(betaw=.03,dxcrtc0=9.e3) parameter(h1=0.33333333) ! progsigma - parameter(dxcrtas=30.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) + parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), & uo(im,km), vo(im,km), qeso(im,km), diff --git a/physics/CONV/SAMF/samfshalcnv.meta b/physics/CONV/SAMF/samfshalcnv.meta index d6fb7960d..c1714801d 100644 --- a/physics/CONV/SAMF/samfshalcnv.meta +++ b/physics/CONV/SAMF/samfshalcnv.meta @@ -504,7 +504,7 @@ [omegain] standard_name = prognostic_updraft_velocity_in_convection long_name = convective updraft velocity - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys @@ -513,7 +513,7 @@ [omegaout] standard_name = updraft_velocity_updated_by_physics long_name = convective updraft velocity updated by physics - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.F90 index 921b45111..632a86597 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.F90 @@ -46,7 +46,7 @@ end subroutine GFS_photochemistry_init !! subroutine GFS_photochemistry_run (dtp, ozphys, oz_phys_2015, oz_phys_2006, con_1ovg, & prsl, dp, ozpl, h2o_phys, h2ophys, h2opl, h2o0, oz0, gt0, do3_dt_prd, do3_dt_ozmx, & - do3_dt_temp, do3_dt_ohoz, errmsg, errflg) + do3_dt_temp, do3_dt_ohoz, dqv_dt_prd, dqv_dt_qvmx, errmsg, errflg) ! Inputs real(kind=kind_phys), intent(in) :: & @@ -73,7 +73,9 @@ subroutine GFS_photochemistry_run (dtp, ozphys, oz_phys_2015, oz_phys_2006, con_ do3_dt_prd, & ! Physics tendency: production and loss effect do3_dt_ozmx, & ! Physics tendency: ozone mixing ratio effect do3_dt_temp, & ! Physics tendency: temperature effect - do3_dt_ohoz ! Physics tendency: overhead ozone effect + do3_dt_ohoz, & ! Physics tendency: overhead ozone effect + dqv_dt_prd, & ! Physics tendency: Climatological net production effect + dqv_dt_qvmx ! Physics tendency: specific humidity effect ! Outputs real(kind=kind_phys), intent(inout), dimension(:,:) :: & @@ -97,7 +99,7 @@ subroutine GFS_photochemistry_run (dtp, ozphys, oz_phys_2015, oz_phys_2006, con_ do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz) endif if (h2o_phys) then - call h2ophys%run(dtp, prsl, h2opl, h2o0) + call h2ophys%run(dtp, prsl, h2opl, h2o0, dqv_dt_prd, dqv_dt_qvmx) endif end subroutine GFS_photochemistry_run diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.meta index f8874fef7..a89773b12 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_photochemistry.meta @@ -193,6 +193,24 @@ kind = kind_phys intent = inout optional = True +[dqv_dt_prd] + standard_name = water_vapor_tendency_due_to_production_and_loss_rate + long_name = water_vapor tendency due to production and loss rate + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True +[dqv_dt_qvmx] + standard_name = water_vapor_tendency_due_to_water_vapor_mixing_ratio + long_name = water_vapor tendency due to water_vapor mixing ratio + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout + optional = True [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90 index 1c78cccb5..06d35c17b 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90 @@ -18,7 +18,8 @@ module GFS_phys_time_vary use module_h2ophys, only: ty_h2ophys use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm, iamin, iamax, jamin, jamax - use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf + use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf, & + read_aerdata_dl, aerinterpol_dl, read_aerdataf_dl use iccn_def, only : ciplin, ccnin, ci_pres use iccninterp, only : read_cidata, setindxci, ciinterpol @@ -78,7 +79,7 @@ end subroutine copy_error !>\section gen_GFS_phys_time_vary_init GFS_phys_time_vary_init General Algorithm !> @{ subroutine GFS_phys_time_vary_init ( & - me, master, ntoz, h2o_phys, iaerclm, iccn, iaermdl, iflip, im, levs, & + me, master, ntoz, h2o_phys, iaerclm, iaermdl, iccn, iflip, im, levs, & nx, ny, idate, xlat_d, xlon_d, & jindx1_o3, jindx2_o3, ddy_o3, jindx1_h, jindx2_h, ddy_h, h2opl,fhour, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & @@ -227,7 +228,12 @@ subroutine GFS_phys_time_vary_init ( ntrcaer = ntrcaerm myerrflg = 0 myerrmsg = 'read_aerdata failed without a message' - call read_aerdata (me,master,iflip,idate,myerrmsg,myerrflg) + if(iaermdl == 1) then + call read_aerdata (me,master,iflip,idate,myerrmsg,myerrflg) + elseif (iaermdl == 6) then + call read_aerdata_dl(me,master,iflip, & + idate,fhour, myerrmsg,myerrflg) + end if call copy_error(myerrmsg, myerrflg, errmsg, errflg) else if(iaermdl ==2 ) then do ix=1,ntrcaerm @@ -351,7 +357,11 @@ subroutine GFS_phys_time_vary_init ( if (iaerclm) then ! This call is outside the OpenMP section, so it should access errmsg & errflg directly. - call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg) + if(iaermdl==1) then + call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg) + elseif (iaermdl==6) then + call read_aerdataf_dl (me, master, iflip, idate, fhour, errmsg, errflg) + end if ! If it is moved to an OpenMP section, it must use myerrmsg, myerrflg, and copy_error. if (errflg/=0) return end if @@ -704,7 +714,7 @@ end subroutine GFS_phys_time_vary_init subroutine GFS_phys_time_vary_timestep_init ( & me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, cplflx, & nsswr, fhswr, lsswr, fhour, & - imfdeepcnv, cal_pre, random_clds, nscyc, ntoz, h2o_phys, iaerclm, iccn, clstp, & + imfdeepcnv, cal_pre, random_clds, nscyc, ntoz, h2o_phys, iaerclm, iaermdl, iccn, clstp, & jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, iflip, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, fn_nml, & @@ -722,7 +732,7 @@ subroutine GFS_phys_time_vary_timestep_init ( ! Interface variables integer, intent(in) :: me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, & - nsswr, imfdeepcnv, iccn, nscyc, ntoz, iflip + nsswr, imfdeepcnv, iccn, nscyc, ntoz, iflip, iaermdl integer, intent(in) :: idate(:) real(kind_phys), intent(in) :: fhswr, fhour logical, intent(in) :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm, cplflx @@ -790,7 +800,7 @@ subroutine GFS_phys_time_vary_timestep_init ( !$OMP parallel num_threads(nthrds) default(none) & !$OMP shared(kdt,nsswr,lsswr,clstp,imfdeepcnv,cal_pre,random_clds) & -!$OMP shared(fhswr,fhour,seed0,cnx,cny,nrcm,wrk,rannie,rndval) & +!$OMP shared(fhswr,fhour,seed0,cnx,cny,nrcm,wrk,rannie,rndval, iaermdl) & !$OMP shared(rann,im,isc,jsc,imap,jmap,ntoz,me,idate,jindx1_o3,jindx2_o3) & !$OMP shared(ozpl,ddy_o3,h2o_phys,jindx1_h,jindx2_h,h2opl,ddy_h,iaerclm,master) & !$OMP shared(levs,prsl,iccn,jindx1_ci,jindx2_ci,ddy_ci,iindx1_ci,iindx2_ci) & @@ -863,23 +873,17 @@ subroutine GFS_phys_time_vary_timestep_init ( rjday = jdoy + jdat(5) / 24. if (rjday < ozphys%time(1)) rjday = rjday + 365. - n2 = ozphys%ntime + 1 - do j=2,ozphys%ntime - if (rjday < ozphys%time(j)) then - n2 = j - exit - endif - enddo - n1 = n2 - 1 - if (n2 > ozphys%ntime) n2 = n2 - ozphys%ntime - !> - Update ozone concentration. if (ntoz > 0) then + call find_photochem_time_index(ozphys%ntime, ozphys%time, rjday, n1, n2) + call ozphys%update_o3prog(jindx1_o3, jindx2_o3, ddy_o3, rjday, n1, n2, ozpl) endif !> - Update stratospheric h2o concentration. if (h2o_phys) then + call find_photochem_time_index(h2ophys%ntime, h2ophys%time, rjday, n1, n2) + call h2ophys%update(jindx1_h, jindx2_h, ddy_h, rjday, n1, n2, h2opl) endif @@ -908,11 +912,19 @@ subroutine GFS_phys_time_vary_timestep_init ( if (iaerclm) then ! aerinterpol is using threading inside, don't ! move into OpenMP parallel section above - call aerinterpol (me, master, nthrds, im, idate, & - fhour, iflip, jindx1_aer, jindx2_aer, & - ddy_aer, iindx1_aer, & - iindx2_aer, ddx_aer, & - levs, prsl, aer_nm, errmsg, errflg) + if (iaermdl==1) then + call aerinterpol (me, master, nthrds, im, idate, & + fhour, iflip, jindx1_aer, jindx2_aer, & + ddy_aer, iindx1_aer, & + iindx2_aer, ddx_aer, & + levs, prsl, aer_nm, errmsg, errflg) + else if (iaermdl==6) then + call aerinterpol_dl (me, master, nthrds, im, idate, & + fhour, iflip, jindx1_aer, jindx2_aer, & + ddy_aer, iindx1_aer, & + iindx2_aer, ddx_aer, & + levs, prsl, aer_nm, errmsg, errflg) + endif if(errflg /= 0) then return endif @@ -933,6 +945,29 @@ subroutine GFS_phys_time_vary_timestep_init ( endif endif + contains + !> Find the time indexes on either side of current time + subroutine find_photochem_time_index(ntime, time, rjday, n1, n2) + implicit none + !> The number of times provided in the parameter file + integer, intent(in) :: ntime + !> The indexes of the parameters just before and after the + !! current time + integer, intent(out) :: n1, n2 + !> The times provided in the parameter file + real, intent(in), dimension(ntime+1) :: time + !> The current time of year + real, intent(in) :: rjday + n2 = ntime + 1 + do j=2,ntime + if (rjday < time(j)) then + n2 = j + exit + endif + enddo + n1 = n2 - 1 + if (n2 > ntime) n2 = n2 - ntime + end subroutine find_photochem_time_index end subroutine GFS_phys_time_vary_timestep_init !> @} diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.meta index 8dba6760e..242cc467d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.meta @@ -1259,6 +1259,13 @@ dimensions = () type = logical intent = in +[iaermdl] + standard_name = control_for_aerosol_radiation_scheme + long_name = control of aerosol scheme in radiation + units = 1 + dimensions = () + type = integer + intent = in [iccn] standard_name = control_for_ice_cloud_condensation_nuclei_forcing long_name = flag for IN and CCN forcing for morrison gettelman microphysics diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 index 35b46618c..d66c1f19f 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 @@ -15,7 +15,8 @@ module GFS_phys_time_vary use module_h2ophys, only: ty_h2ophys use aerclm_def, only : aerin, aer_pres, ntrcaer, ntrcaerm, iamin, iamax, jamin, jamax - use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf + use aerinterp, only : read_aerdata, setindxaer, aerinterpol, read_aerdataf, & + read_aerdata_dl, aerinterpol_dl, read_aerdataf_dl use iccn_def, only : ciplin, ccnin, ci_pres use iccninterp, only : read_cidata, setindxci, ciinterpol @@ -58,7 +59,7 @@ module GFS_phys_time_vary !>\section gen_GFS_phys_time_vary_init GFS_phys_time_vary_init General Algorithm !! @{ subroutine GFS_phys_time_vary_init ( & - me, master, ntoz, h2o_phys, iaerclm, iccn, iflip, im, nx, ny, idate, xlat_d, xlon_d, & + me, master, ntoz, h2o_phys, iaerclm,iaermdl, iccn, iflip, im, nx, ny, idate, xlat_d, xlon_d, & jindx1_o3, jindx2_o3, ddy_o3, ozphys, h2ophys, jindx1_h, jindx2_h, ddy_h, h2opl,fhour, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap, & @@ -77,7 +78,7 @@ subroutine GFS_phys_time_vary_init ( implicit none ! Interface variables - integer, intent(in) :: me, master, ntoz, iccn, iflip, im, nx, ny + integer, intent(in) :: me, master, ntoz, iccn, iflip, im, nx, ny, iaermdl logical, intent(in) :: h2o_phys, iaerclm, lsm_cold_start integer, intent(in) :: idate(:), iopt_lake, iopt_lake_clm, iopt_lake_flake real(kind_phys), intent(in) :: fhour, lakefrac_threshold, lakedepth_threshold @@ -208,7 +209,11 @@ subroutine GFS_phys_time_vary_init ( ! If iaerclm is .true., then ntrcaer == ntrcaerm ntrcaer = size(aer_nm, dim=3) ! Read aerosol climatology - call read_aerdata (me,master,iflip,idate,errmsg,errflg) + if(iaermdl==1) then + call read_aerdata (me,master,iflip,idate,errmsg,errflg) + elseif(iaermdl==6) then + call read_aerdata_dl (me,master,iflip,idate,fhour,errmsg,errflg) + end if endif if (errflg /= 0) return else @@ -312,9 +317,12 @@ subroutine GFS_phys_time_vary_init ( endif if (errflg/=0) return - if (iaerclm) then - call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg) + if (iaermdl==1) then + call read_aerdataf (me, master, iflip, idate, fhour, errmsg, errflg) + elseif (iaermdl==6) then + call read_aerdataf_dl (me, master, iflip, idate, fhour, errmsg, errflg) + end if if (errflg/=0) return end if @@ -646,7 +654,7 @@ end subroutine GFS_phys_time_vary_init !! @{ subroutine GFS_phys_time_vary_timestep_init ( & me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, nsswr, fhswr, lsswr, fhour, & - imfdeepcnv, cal_pre, random_clds, ozphys, h2ophys, ntoz, h2o_phys, iaerclm, iccn, clstp, & + imfdeepcnv, cal_pre, random_clds, ozphys, h2ophys, ntoz, h2o_phys, iaerclm, iaermdl, iccn, clstp, & jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, iflip, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, & @@ -657,7 +665,7 @@ subroutine GFS_phys_time_vary_timestep_init ( ! Interface variables integer, intent(in) :: me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, & - nsswr, imfdeepcnv, iccn, ntoz, iflip + nsswr, imfdeepcnv, iccn, ntoz, iflip, iaermdl integer, intent(in) :: idate(:) real(kind_phys), intent(in) :: fhswr, fhour logical, intent(in) :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm @@ -769,23 +777,17 @@ subroutine GFS_phys_time_vary_timestep_init ( rjday = jdoy + jdat(5) / 24. if (rjday < ozphys%time(1)) rjday = rjday + 365. - n2 = ozphys%ntime + 1 - do j=2,ozphys%ntime - if (rjday < ozphys%time(j)) then - n2 = j - exit - endif - enddo - n1 = n2 - 1 - if (n2 > ozphys%ntime) n2 = n2 - ozphys%ntime - !> - Update ozone concentration. if (ntoz > 0) then + call find_photochemistry_index(ozphys%ntime, ozphys%time, rjday, n1, n2) + call ozphys%update_o3prog(jindx1_o3, jindx2_o3, ddy_o3, rjday, n1, n2, ozpl) endif !> - Update stratospheric h2o concentration. if (h2o_phys) then + call find_photochemistry_index(h2ophys%ntime, h2ophys%time, rjday, n1, n2) + call h2ophys%update(jindx1_h, jindx2_h, ddy_h, rjday, n1, n2, h2opl) endif @@ -809,13 +811,18 @@ subroutine GFS_phys_time_vary_timestep_init ( if (iaerclm) then ! aerinterpol is using threading inside, don't ! move into OpenMP parallel section above - call aerinterpol (me, master, nthrds, im, idate, & - fhour, iflip, jindx1_aer, jindx2_aer, & - ddy_aer, iindx1_aer, & - iindx2_aer, ddx_aer, & - levs, prsl, aer_nm, errmsg, errflg) - if(errflg /= 0) then - return + if (iaermdl==1) then + call aerinterpol (me, master, nthrds, im, idate, & + fhour, iflip, jindx1_aer, jindx2_aer, & + ddy_aer, iindx1_aer, & + iindx2_aer, ddx_aer, & + levs, prsl, aer_nm, errmsg, errflg) + elseif (iaermdl==6) then + call aerinterpol_dl (me, master, nthrds, im, idate, & + fhour, iflip, jindx1_aer, jindx2_aer, & + ddy_aer, iindx1_aer, & + iindx2_aer, ddx_aer, & + levs, prsl, aer_nm, errmsg, errflg) endif endif @@ -834,6 +841,29 @@ subroutine GFS_phys_time_vary_timestep_init ( ! endif ! endif + contains + !> Find the time indexes on either side of current time + subroutine find_photochemistry_index(ntime, time, rjday, n1, n2) + implicit none + !> The number of times provided in the parameter file + integer, intent(in) :: ntime + !> The indexes of the parameters just before and after the + !! current time + integer, intent(out) :: n1, n2 + !> The times provided in the parameter file + real, intent(in), dimension(ntime+1) :: time + !> The current time of year + real, intent(in) :: rjday + n2 = ntime + 1 + do j=2,ntime + if (rjday < time(j)) then + n2 = j + exit + endif + enddo + n1 = n2 - 1 + if (n2 > ntime) n2 = n2 - ntime + end subroutine find_photochemistry_index end subroutine GFS_phys_time_vary_timestep_init !! @} diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta index 5dae234ea..9bc9c39b9 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta @@ -50,6 +50,13 @@ dimensions = () type = logical intent = in +[iaermdl] + standard_name = control_for_aerosol_radiation_scheme + long_name = control of aerosol scheme in radiation + units = 1 + dimensions = () + type = integer + intent = in [iccn] standard_name = control_for_ice_cloud_condensation_nuclei_forcing long_name = flag for IN and CCN forcing for morrison gettelman microphysics @@ -1216,6 +1223,13 @@ dimensions = () type = logical intent = in +[iaermdl] + standard_name = control_for_aerosol_radiation_scheme + long_name = control of aerosol scheme in radiation + units = 1 + dimensions = () + type = integer + intent = in [iccn] standard_name = control_for_ice_cloud_condensation_nuclei_forcing long_name = flag for IN and CCN forcing for morrison gettelman microphysics diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.F90 index 37ac9b320..7bb3c7e6d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.F90 @@ -18,6 +18,7 @@ module GFS_physics_post subroutine GFS_physics_post_run(nCol, nLev, ntoz, ntracp100, nprocess, nprocess_summed, & dtidx, is_photochem, ldiag3d, ip_physics, ip_photochem, ip_prod_loss, ip_ozmix, & ip_temp, ip_overhead_ozone, do3_dt_prd, do3_dt_ozmx, do3_dt_temp, do3_dt_ohoz, & + ntqv, dqv_dt_prd, dqv_dt_qvmx, & dtend, errmsg, errflg) ! Inputs @@ -25,6 +26,7 @@ subroutine GFS_physics_post_run(nCol, nLev, ntoz, ntracp100, nprocess, nprocess_ nCol, & !< Horizontal dimension nLev, & !< Number of vertical layers ntoz, & !< Index for ozone mixing ratio + ntqv, & !< Index for water vapor mixing ratio ntracp100, & !< Number of tracers plus 100 nprocess, & !< Number of processes that cause changes in state variables nprocess_summed,& !< Number of causes in dtidx per tracer summed for total physics tendency @@ -46,7 +48,9 @@ subroutine GFS_physics_post_run(nCol, nLev, ntoz, ntracp100, nprocess, nprocess_ do3_dt_prd, & !< Physics tendency: production and loss effect do3_dt_ozmx, & !< Physics tendency: ozone mixing ratio effect do3_dt_temp, & !< Physics tendency: temperature effect - do3_dt_ohoz !< Physics tendency: overhead ozone effect + do3_dt_ohoz, & !< Physics tendency: overhead ozone effect + dqv_dt_prd, & !< Physics tendency: climatological net production effect + dqv_dt_qvmx !< Physics tendency: water vapor mixing ratio effect ! Outputs real(kind=kind_phys), intent(inout), dimension(:,:,:), optional :: & @@ -104,6 +108,32 @@ subroutine GFS_physics_post_run(nCol, nLev, ntoz, ntracp100, nprocess, nprocess_ call sum_it(ichem, itrac, is_photochem) endif + ! ####################################################################################### + ! + ! Water vapor photochemistry diagnostics + ! + ! ####################################################################################### + idtend = dtidx(100+ntqv, ip_prod_loss) + if (idtend >= 1 .and. associated(dqv_dt_prd)) then + dtend(:, :, idtend) = dtend(:, :, idtend) + dqv_dt_prd + endif + ! + idtend = dtidx(100+ntqv,ip_ozmix) + if (idtend >= 1 .and. associated(dqv_dt_qvmx)) then + dtend(:, :, idtend) = dtend(:, :, idtend) + dqv_dt_qvmx + endif + + ! ####################################################################################### + ! + ! Total photochemical tendencies + ! + ! ####################################################################################### + itrac = ntqv + 100 + ichem = dtidx(itrac, ip_photochem) + if (ichem >= 1) then + call sum_it(ichem, itrac, is_photochem) + endif + ! ####################################################################################### ! ! Total (physics) tendencies diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.meta index c42da0166..823870a0e 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_physics_post.meta @@ -157,6 +157,31 @@ kind = kind_phys intent = in optional = True +[ntqv] + standard_name = index_of_specific_humidity_in_tracer_concentration_array + long_name = tracer index for water_vapor mixing ratio + units = index + dimensions = () + type = integer + intent = in +[dqv_dt_prd] + standard_name = water_vapor_tendency_due_to_production_and_loss_rate + long_name = water_vapor tendency due to production and loss rate + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = True +[dqv_dt_qvmx] + standard_name = water_vapor_tendency_due_to_water_vapor_mixing_ratio + long_name = water_vapor tendency due to water_vapor mixing ratio + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in + optional = True [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.F90 index cc8e950e8..804638f8a 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_setup.F90 @@ -184,12 +184,6 @@ subroutine GFS_rrtmg_setup_init ( si, levr, ictm, isol, solar_file, ico2, & iaerflg = mod(iaer, 1000) endif iaermdl = iaer/1000 ! control flag for aerosol scheme selection - if ( iaermdl < 0 .or. (iaermdl>2 .and. iaermdl/=5) ) then - print *, ' Error -- IAER flag is incorrect, Abort' - errflg = 1 - errmsg = 'ERROR(GFS_rrtmg_setup): IAER flag is incorrect' - return - endif ! --- assign initial permutation seed for mcica cloud-radiation if ( isubcsw>0 .or. isubclw>0 ) then diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 index 0ed936410..7c95d726d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 @@ -88,11 +88,6 @@ subroutine GFS_rrtmgp_setup_init(do_RRTMGP, imp_physics, imp_physics_fer_hires, iaerflg = mod(iaer, 1000) endif iaermdl = iaer/1000 ! control flag for aerosol scheme selection - if ( iaermdl < 0 .or. (iaermdl>2 .and. iaermdl/=5) ) then - errmsg = trim(errmsg) // ' Error -- IAER flag is incorrect, Abort' - errflg = 1 - return - endif ! Assign initial permutation seed for mcica cloud-radiation if ( isubc_sw>0 .or. isubc_lw>0 ) then diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta index 2f267aa75..4cf339e4d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta @@ -273,7 +273,7 @@ [omegain] standard_name = prognostic_updraft_velocity_in_convection long_name = convective updraft velocity - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys @@ -282,7 +282,7 @@ [omegaout] standard_name = updraft_velocity_updated_by_physics long_name = convective updraft velocity updated by physics - units = frac + units = Pa s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys diff --git a/physics/MP/Morrison_Gettelman/aerclm_def.F b/physics/MP/Morrison_Gettelman/aerclm_def.F index b6760f30c..3bfa6791a 100644 --- a/physics/MP/Morrison_Gettelman/aerclm_def.F +++ b/physics/MP/Morrison_Gettelman/aerclm_def.F @@ -3,15 +3,17 @@ module aerclm_def implicit none integer, parameter :: levsaer=72, ntrcaerm=15, timeaer=2 - integer :: latsaer, lonsaer, ntrcaer, levsw - integer :: n1sv, n2sv + integer :: latsaer, lonsaer, ntrcaer, levsw, tsaer + integer :: n1sv, n2sv, t1sv, t2sv integer :: iamin, iamax, jamin, jamax character*10 :: specname(ntrcaerm) + character*50 :: fname_dl real (kind=kind_phys):: aer_time(13) real (kind=kind_phys), allocatable, dimension(:) :: aer_lat real (kind=kind_phys), allocatable, dimension(:) :: aer_lon + real (kind=kind_phys), allocatable, dimension(:) :: aer_t real (kind=kind_io4), allocatable, dimension(:,:,:,:) :: aer_pres real (kind=kind_io4), allocatable, dimension(:,:,:,:,:) :: aerin diff --git a/physics/MP/Morrison_Gettelman/aerinterp.F90 b/physics/MP/Morrison_Gettelman/aerinterp.F90 index 97113a4e1..15963db7d 100644 --- a/physics/MP/Morrison_Gettelman/aerinterp.F90 +++ b/physics/MP/Morrison_Gettelman/aerinterp.F90 @@ -9,16 +9,17 @@ module aerinterp implicit none - private read_netfaer + private read_netfaer, read_netfaer_dl, fdnx_fname public :: read_aerdata, setindxaer, aerinterpol,read_aerdataf + public :: read_aerdata_dl, aerinterpol_dl,read_aerdataf_dl contains logical function netcdf_check(status, errmsg, errflg, why) use netcdf implicit none - character(len=*), intent(inout) :: errmsg + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg integer, intent(in) :: status character(len=*), intent(in) :: why @@ -34,21 +35,407 @@ logical function netcdf_check(status, errmsg, errflg, why) endif END function netcdf_check - - SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) - use machine, only: kind_phys, kind_io4, kind_io8 +!!!!!!! + SUBROUTINE read_aerdata_dl (me, master, iflip, idate, FHOUR, errmsg, errflg) + use machine, only: kind_phys, kind_io4, kind_dbl_prec use aerclm_def use netcdf +!--- in/out + integer, intent(in) :: me, master, iflip, idate(4) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + real(kind=kind_phys), intent(in) :: fhour +!--- locals + integer :: ncid, varid, ndims, hmx + integer :: i, j, k, n, ii, imon, klev + character :: fname*50, mn*2, vname*10, dy*2, myr*4 + logical :: file_exist + integer :: dimids(NF90_MAX_VAR_DIMS) + integer :: dimlen(NF90_MAX_VAR_DIMS) + integer IDAT(8),JDAT(8) + real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) + integer jdow, jdoy, jday + + errflg = 0 + errmsg = ' ' + +! +!! =================================================================== + if (me == master) then + if ( iflip == 0 ) then ! data from toa to sfc + print *, "GFS is top-down" + else + print *, "GFS is bottom-up" + endif + endif +!! found first day needed to interpolated + IDAT = 0 + IDAT(1) = IDATE(4) + IDAT(2) = IDATE(2) + IDAT(3) = IDATE(3) + IDAT(5) = IDATE(1) + RINC = 0. + RINC(2) = FHOUR + CALL W3MOVDAT(RINC,IDAT,JDAT) +! + jdow = 0 + jdoy = 0 + jday = 0 + call w3doxdat(jdat,jdow,jdoy,jday) + +! +!! =================================================================== +!! check if all necessary files exist +!! =================================================================== + write(myr,'(i4.4)') jdat(1) + write(mn,'(i2.2)') jdat(2) + write(dy,'(i2.2)') jdat(3) + fname=trim("merra2_"//myr//mn//dy//".nc") + inquire (file = fname, exist = file_exist) + if (.not. file_exist) then + errmsg = 'Error in read_aerdata: file ' // trim(fname) // ' not found' + errflg = 1 + return + endif +! +!! =================================================================== +!! fetch dim spec and lat/lon from m01 data set +!! =================================================================== + ncid = -1 + if(.not.netcdf_check(nf90_open(fname , nf90_NOWRITE, ncid), & + errmsg, errflg, 'open '//trim(fname))) then + return + endif + + vname = trim(specname(1)) + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), & + errmsg, errflg, 'find id of '//trim(vname)//' var')) then + return + endif + ndims = 0 + if(.not.netcdf_check(nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dimids), & + errmsg, errflg, 'inquire details about '//trim(vname)//' var')) then + return + endif + do i=1,ndims + if(.not.netcdf_check(nf90_inquire_dimension(ncid, dimids(i), len=dimlen(i)), & + errmsg, errflg, 'inquire details about dimension')) then + return + endif + enddo + +! specify latsaer, lonsaer, hmx + lonsaer = dimlen(1) + latsaer = dimlen(2) + levsw = dimlen(3) + tsaer = dimlen(4) + + if(me==master) then + print *, 'MERRA2 dim: ',dimlen(1:ndims) + endif + +! allocate arrays + + if (.not. allocated(aer_lat)) then + allocate(aer_lat(latsaer)) + allocate(aer_lon(lonsaer)) + allocate(aer_t(tsaer)) + endif + +! construct lat/lon array + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, 'lat', varid), & + errmsg, errflg, 'find id of lat var')) then + return + endif + aer_lat = 0 + if(.not.netcdf_check(nf90_get_var(ncid, varid, aer_lat, (/ 1, 1, 1 /), (/latsaer, 1, 1/)), & + errmsg, errflg, 'read lat var')) then + return + endif + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, 'lon', varid), & + errmsg, errflg, 'find id of lon var')) then + return + endif + aer_lon = 0 + if(.not.netcdf_check(nf90_get_var(ncid, varid, aer_lon, (/ 1, 1, 1 /), (/lonsaer, 1, 1/)), & + errmsg, errflg, 'read lon var')) then + return + endif + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, 'time', varid), & + errmsg, errflg, 'find id of time var')) then + return + endif + aer_t = 0 + if(.not.netcdf_check(nf90_get_var(ncid, varid, aer_t, (/ 1, 1, 1 /), (/tsaer, 1, 1/)), & + errmsg, errflg, 'read t var')) then + return + endif + if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg, 'close '//trim(fname))) then + return + endif + END SUBROUTINE read_aerdata_dl +! +!********************************************************************** + SUBROUTINE read_aerdataf_dl ( me, master, iflip, idate, FHOUR, errmsg, errflg) + use machine, only: kind_phys, kind_dbl_prec + use aerclm_def + !--- in/out integer, intent(in) :: me, master, iflip, idate(4) character(len=*), intent(inout) :: errmsg integer, intent(inout) :: errflg + real(kind=kind_phys), intent(in) :: fhour +!--- locals + integer :: i, j, k, n, ii, imon, klev, n1, n2 + logical :: file_exist, fd_upb + integer IDAT(8),JDAT(8) + real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) + integer jdow, jdoy, jday + character myr*4, mn*2, dy*2 + + integer, allocatable :: invardims(:) +! + if (.not. allocated(aerin)) then + allocate(aerin(iamin:iamax,jamin:jamax,levsaer,ntrcaerm,timeaer)) + allocate(aer_pres(iamin:iamax,jamin:jamax,levsaer,timeaer)) + endif + +! allocate local working arrays +!! found interpolation months + IDAT = 0 + IDAT(1) = IDATE(4) + IDAT(2) = IDATE(2) + IDAT(3) = IDATE(3) + IDAT(5) = IDATE(1) + RINC = 0. + RINC(2) = FHOUR + CALL W3MOVDAT(RINC,IDAT,JDAT) +! + jdow = 0 + jdoy = 0 + jday = 0 + call w3doxdat(jdat,jdow,jdoy,jday) + write(myr,'(i4.4)') jdat(1) + write(mn,'(i2.2)') jdat(2) + write(dy,'(i2.2)') jdat(3) + fname_dl="merra2_"//myr//mn//dy//".nc" +! rjday is the minutes in a day + rjday = jdat(5)*60+jdat(6)+jdat(7)/60. +! +! n1sv saves the jdat(3), n2sv is the up boundary index + n1sv=jdat(3) + fd_upb=.false. + do j=2, tsaer + if ( aer_t(j)> rjday) then + n2sv = j + t2sv = aer_t(j) + fd_upb=.true. + exit + endif + enddo + if(fd_upb) then + t1sv = aer_t(j-1) + call read_netfaer_dl(fname_dl, j-1, iflip, 1, errmsg, errflg) + call read_netfaer_dl(fname_dl, n2sv, iflip, 2, errmsg, errflg) + else + t1sv = aer_t(tsaer) + call read_netfaer_dl(fname_dl, tsaer, iflip, 1, errmsg, errflg) + n2sv=1 + t2sv=1440. + call fdnx_fname (jdat(1), jdat(2),jdat(3),fname_dl) + call read_netfaer_dl(fname_dl, n2sv, iflip, 2, errmsg, errflg) + end if + END SUBROUTINE read_aerdataf_dl +!********************************************************************** +! + SUBROUTINE aerinterpol_dl( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, & + ddy,iindx1,iindx2,ddx,lev,prsl,aerout, errmsg,errflg) +! + use machine, only: kind_phys, kind_dbl_prec + use aerclm_def + + implicit none + integer, intent(out) :: errflg + character(*), intent(out) :: errmsg + integer, intent(in) :: iflip + integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii, klev + real(kind=kind_phys) fhour,temj, tx1, tx2,temi, tem, tem1, tem2 + real(kind=kind_phys), dimension(npts) :: temij,temiy,temjx,ddxy + +! + + integer JINDX1(npts), JINDX2(npts), iINDX1(npts), iINDX2(npts) + integer me,idate(4), master, nthrds + integer IDAT(8),JDAT(8) +! + real(kind=kind_phys) DDY(npts), ddx(npts),ttt + real(kind=kind_phys) aerout(npts,lev,ntrcaer) + real(kind=kind_phys) aerpm(npts,levsaer,ntrcaer) + real(kind=kind_phys) prsl(npts,lev), aerpres(npts,levsaer) + real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) + integer jdow, jdoy, jday + character myr*4, mn*2, dy*2 +! + errflg = 0 + errmsg = ' ' + IDAT = 0 + IDAT(1) = IDATE(4) + IDAT(2) = IDATE(2) + IDAT(3) = IDATE(3) + IDAT(5) = IDATE(1) + RINC = 0. + RINC(2) = FHOUR + CALL W3MOVDAT(RINC,IDAT,JDAT) +! + jdow = 0 + jdoy = 0 + jday = 0 + call w3doxdat(jdat,jdow,jdoy,jday) +! rjday is the minutes in a day + rjday = jdat(5)*60+jdat(6)+jdat(7)/60. + if(rjday >= t2sv .or. jdat(3).ne.n1sv) then !!need to either to read in a record or open a new file + call read_netfaer_dl(fname_dl,n2sv, iflip, 1, errmsg, errflg) + end if +!! =================================================================== + if(jdat(3).ne.n1sv) then ! a new day is produced from n2sv=1440 + n1sv=jdat(3) + t1sv=aer_t(1) + n2sv=2 + t2sv=aer_t(n2sv) + write(myr,'(i4.4)') jdat(1) + write(mn,'(i2.2)') jdat(2) + write(dy,'(i2.2)') jdat(3) + fname_dl="merra2_"//myr//mn//dy//".nc" + call read_netfaer_dl(fname_dl,n2sv, iflip, 2, errmsg, errflg) + else if (rjday >= t2sv) then + if(t2sv < aer_t(tsaer)) then + n1sv=jdat(3) + t1sv=t2sv + n2sv=n2sv+1 + t2sv=aer_t(n2sv) + write(myr,'(i4.4)') jdat(1) + write(mn,'(i2.2)') jdat(2) + write(dy,'(i2.2)') jdat(3) + fname_dl="merra2_"//myr//mn//dy//".nc" + call read_netfaer_dl(fname_dl,n2sv, iflip, 2, errmsg, errflg) + else !! need to read a new file + n1sv=jdat(3) + t1sv=aer_t(tsaer) + n2sv=1 + t2sv=1440. + call fdnx_fname (jdat(1), jdat(2),jdat(3),fname_dl) + call read_netfaer_dl(fname_dl, n2sv, iflip, 2, errmsg, errflg) + end if + end if +! + tx1 = (t2sv - rjday) / (t2sv - t1sv) + tx2 = 1.0 - tx1 + if (n2 > 12) n2 = n2 -12 + + do j=1,npts + TEMJ = 1.0 - DDY(J) + TEMI = 1.0 - DDX(J) + temij(j) = TEMI*TEMJ + temiy(j) = TEMI*DDY(j) + temjx(j) = TEMJ*DDX(j) + ddxy(j) = DDX(j)*DDY(J) + enddo + +#ifndef __GFORTRAN__ +!$OMP parallel num_threads(nthrds) default(none) & +!$OMP shared(npts,ntrcaer,aerin,aer_pres,prsl) & +!$OMP shared(ddx,ddy,jindx1,jindx2,iindx1,iindx2) & +!$OMP shared(aerpm,aerpres,aerout,lev,nthrds) & +!$OMP shared(temij,temiy,temjx,ddxy,tx1,tx2) & +!$OMP private(l,j,k,ii,i1,i2,j1,j2,tem,tem1,tem2) + +!$OMP do +#endif + DO L=1,levsaer + DO J=1,npts + J1 = JINDX1(J) + J2 = JINDX2(J) + I1 = IINDX1(J) + I2 = IINDX2(J) + DO ii=1,ntrcaer + aerpm(j,L,ii) = & + tx1*(TEMIJ(j)*aerin(I1,J1,L,ii,1)+DDXY(j)*aerin(I2,J2,L,ii,1) & + +TEMIY(j)*aerin(I1,J2,L,ii,1)+temjx(j)*aerin(I2,J1,L,ii,1))& + +tx2*(TEMIJ(j)*aerin(I1,J1,L,ii,2)+DDXY(j)*aerin(I2,J2,L,ii,2) & + +TEMIY(j)*aerin(I1,J2,L,ii,2)+temjx(j)*aerin(I2,J1,L,ii,2)) + ENDDO + + aerpres(j,L) = & + tx1*(TEMIJ(j)*aer_pres(I1,J1,L,1)+DDXY(j)*aer_pres(I2,J2,L,1) & + +TEMIY(j)*aer_pres(I1,J2,L,1)+temjx(j)*aer_pres(I2,J1,L,1))& + +tx2*(TEMIJ(j)*aer_pres(I1,J1,L,2)+DDXY(j)*aer_pres(I2,J2,L,2) & + +TEMIY(j)*aer_pres(I1,J2,L,2)+temjx(j)*aer_pres(I2,J1,L,2)) + ENDDO + ENDDO +#ifndef __GFORTRAN__ +!$OMP end do + +! don't flip, input is the same direction as GFS (bottom-up) +!$OMP do +#endif + DO J=1,npts + DO L=1,lev + if(prsl(j,L) >= aerpres(j,1)) then + DO ii=1, ntrcaer + aerout(j,L,ii) = aerpm(j,1,ii) !! sfc level + ENDDO + else if(prsl(j,L) <= aerpres(j,levsaer)) then + DO ii=1, ntrcaer + aerout(j,L,ii) = aerpm(j,levsaer,ii) !! toa top + ENDDO + else + DO k=1, levsaer-1 !! from sfc to toa + IF(prsl(j,L) <= aerpres(j,k) .and. prsl(j,L)>aerpres(j,k+1)) then + i1 = k + i2 = min(k+1,levsaer) + exit + ENDIF + ENDDO + tem = 1.0 / (aerpres(j,i1) - aerpres(j,i2)) + tem1 = (prsl(j,L) - aerpres(j,i2)) * tem + tem2 = (aerpres(j,i1) - prsl(j,L)) * tem + DO ii = 1, ntrcaer + aerout(j,L,ii) = aerpm(j,i1,ii)*tem1 + aerpm(j,i2,ii)*tem2 + ENDDO + endif + ENDDO !L-loop + ENDDO !J-loop +#ifndef __GFORTRAN__ +!$OMP end do + +!$OMP end parallel +#endif + + RETURN + END SUBROUTINE aerinterpol_dl + + SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) + use machine, only: kind_phys, kind_io4 + use aerclm_def + use netcdf + +!--- in/out + integer, intent(in) :: me, master, iflip, idate(4) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg !--- locals integer :: ncid, varid, ndims, hmx integer :: i, j, k, n, ii, imon, klev - character :: fname*50, mn*2, vname*10 + character :: fname*50, myr*4, mn*2, dy*2,vname*10 logical :: file_exist integer :: dimids(NF90_MAX_VAR_DIMS) integer :: dimlen(NF90_MAX_VAR_DIMS) @@ -67,7 +454,7 @@ SUBROUTINE read_aerdata (me, master, iflip, idate, errmsg, errflg) endif ! !! =================================================================== -!! check if all necessary files exist +!! check if one file exist !! =================================================================== do imon = 1, 12 write(mn,'(i2.2)') imon @@ -278,8 +665,8 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, use aerclm_def implicit none - integer, intent(inout) :: errflg - character(*), intent(inout) :: errmsg + integer, intent(out) :: errflg + character(*), intent(out) :: errmsg integer, intent(in) :: iflip integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i,ii, klev real(kind=kind_phys) fhour,temj, tx1, tx2,temi, tem, tem1, tem2 @@ -426,14 +813,162 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, RETURN END SUBROUTINE aerinterpol +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! private subroutines + subroutine fdnx_fname(lyear, lmn, ldy, fname) + integer, intent(inout) ::lyear, lmn, ldy + character, intent(out) :: fname*50 + integer, dimension(12) :: mndy + character myr*4, mn*2, dy*2 + data mndy/31, 28, 31,30,31,30,31,31,30,31,30,31/ + ldy=ldy+1 + if(lmn==12) then + if(ldy>mndy(12)) then + ldy=1 + lmn=1 + lyear=lyear+1 + end if + else if(lmn==2) then + if (mod(lyear,4)==0) then + if(ldy>mndy(2)+1) then + ldy=1 + lmn=3 + end if + else + if(ldy>mndy(2)) then + ldy=1 + lmn=3 + end if + end if + else + if(ldy>mndy(lmn)) then + ldy=1 + lmn=lmn+1 + end if + end if + write(myr,'(i4.4)') lyear + write(mn,'(i2.2)') lmn + write(dy,'(i2.2)') ldy + fname="merra2_"//myr//mn//dy//".nc" + RETURN + END SUBROUTINE fdnx_fname + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine read_netfaer_dl(fname, nf, iflip,nt, errmsg,errflg) + use machine, only: kind_phys, kind_io4 + use aerclm_def + use netcdf + integer, intent(in) :: iflip, nf, nt + character,intent(in) :: fname*50 + integer, intent(out) :: errflg + character(*), intent(out) :: errmsg + integer :: ncid, varid, i,j,k,ii,klev + character :: vname*10 + real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff + real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx + real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp + integer lstart(4), lcount(4) +!! =================================================================== + allocate (buff(lonsaer, latsaer, levsw)) + allocate (pres_tmp(lonsaer, levsw)) + allocate (buffx(lonsaer, latsaer, levsw, 1)) + + errflg = 0 + errmsg = ' ' + buff = 0 + pres_tmp = 0 + buffx = 0 + + ncid = -1 + if(.not.netcdf_check(nf90_open(trim(fname), nf90_NOWRITE, ncid), & + errmsg, errflg, 'open '//trim(fname))) then + return + endif + lstart=(/1,1,1,nf/) + lcount=(/lonsaer, latsaer, levsw,1/) +! ====> construct 3-d pressure array (Pa) + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, "DELP", varid), & + errmsg, errflg, 'find id of DELP var')) then + return + endif + if(.not.netcdf_check(nf90_get_var(ncid, varid, buff, lstart, lcount), & + errmsg, errflg, 'read DELP var')) then + return + endif + + do j = jamin, jamax + do i = iamin, iamax +! constract pres_tmp (top-down), note input is top-down + pres_tmp(i,1) = 0. + do k=2, levsw + pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) + enddo !k-loop + enddo !i-loop (lon) + +! extract pres_tmp to fill aer_pres (in Pa) + do k = 1, levsaer + if ( iflip == 0 ) then ! data from toa to sfc + klev = k + else ! data from sfc to top + klev = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aer_pres(i,j,k,nt) = 1.d0*pres_tmp(i,klev) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + +! ====> construct 4-d aerosol array (kg/kg) +! merra2 data is top down +! for GFS, iflip 0: toa to sfc; 1: sfc to toa + DO ii = 1, ntrcaerm + vname=trim(specname(ii)) + varid = -1 + if(.not.netcdf_check(nf90_inq_varid(ncid, vname, varid), & + errmsg, errflg, 'get id of '//trim(vname)//' var')) then + return + endif + if(.not.netcdf_check(nf90_get_var(ncid, varid, buffx, lstart, lcount), & + errmsg, errflg, 'read '//trim(vname)//' var')) then + return + endif + + do j = jamin, jamax + do k = 1, levsaer +! input is from toa to sfc + if ( iflip == 0 ) then ! data from toa to sfc + klev = k + else ! data from sfc to top + klev = ( levsw - k ) + 1 + endif + do i = iamin, iamax + aerin(i,j,k,ii,nt) = 1.d0*buffx(i,j,klev,1) + if(aerin(i,j,k,ii,nt) < 0 .or. aerin(i,j,k,ii,nt) > 1.) then + aerin(i,j,k,ii,nt) = 1.e-15 + endif + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + + ENDDO ! ii-loop (ntracaerm) + +! close the file + if(.not.netcdf_check(nf90_close(ncid), errmsg, errflg, 'close '//trim(fname))) then + return + endif + deallocate (buff, pres_tmp) + deallocate (buffx) + END SUBROUTINE read_netfaer_dl +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_netfaer(nf, iflip,nt, errmsg,errflg) - use machine, only: kind_phys, kind_io4, kind_io8 + use machine, only: kind_phys, kind_io4 use aerclm_def use netcdf integer, intent(in) :: iflip, nf, nt - integer, intent(inout) :: errflg - character(*), intent(inout) :: errmsg + integer, intent(out) :: errflg + character(*), intent(out) :: errmsg integer :: ncid, varid, i,j,k,ii,klev character :: fname*50, mn*2, vname*10 real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff diff --git a/physics/Radiation/radiation_aerosols.f b/physics/Radiation/radiation_aerosols.f index cd5dd49f1..af39ae212 100644 --- a/physics/Radiation/radiation_aerosols.f +++ b/physics/Radiation/radiation_aerosols.f @@ -637,7 +637,7 @@ subroutine aer_init & ! --- outputs: & errflg, errmsg) - elseif ( iaermdl==1 .or. iaermdl==2 ) then ! gocart clim/prog scheme + elseif ( iaermdl==1 .or. iaermdl==2 .or. iaermdl==6) then ! gocart clim/prog scheme call gocart_aerinit & ! --- inputs: @@ -726,7 +726,10 @@ subroutine wrt_aerlog(iaermdl, iaerflg, lalw1bd, errflg, errmsg) print *,' - Using OPAC-seasonal climatology for tropospheric', & & ' aerosol effect' elseif ( iaermdl == 1 ) then - print *,' - Using GOCART-climatology for tropospheric', & + print *,' - Using MERRA2-climatology for tropospheric', & + & ' aerosol effect' + elseif ( iaermdl == 6 ) then + print *,' - Using MERRA2 3 hourly aerosol for tropospheric', & & ' aerosol effect' elseif ( iaermdl == 2 ) then print *,' - Using GOCART-prognostic aerosols for tropospheric', & @@ -2405,7 +2408,7 @@ subroutine setaer & & ) ! - elseif ( iaermdl==1 .or. iaermdl==2) then ! use gocart aerosols + elseif ( iaermdl==1 .or. iaermdl==2 .or. iaermdl==6) then ! use gocart aerosols call aer_property_gocart & ! --- inputs: diff --git a/physics/photochem/module_h2ophys.F90 b/physics/photochem/module_h2ophys.F90 index 61a01d7b4..953232eb3 100644 --- a/physics/photochem/module_h2ophys.F90 +++ b/physics/photochem/module_h2ophys.F90 @@ -149,7 +149,7 @@ end subroutine update ! ######################################################################################### ! Procedure (type-bound) for NRL stratospheric h2o photochemistry physics. ! ######################################################################################### - subroutine run(this, dt, p, h2opltc, h2o) + subroutine run(this, dt, p, h2opltc, h2o, dqv_dt_prd, dqv_dt_qv) class(ty_h2ophys), intent(in) :: this real(kind_phys), intent(in) :: & dt ! Model timestep (sec) @@ -159,6 +159,9 @@ subroutine run(this, dt, p, h2opltc, h2o) h2opltc ! h2o forcing data real(kind_phys), intent(inout), dimension(:,:) :: & h2o ! h2o concentration (updated) + real(kind_phys), intent(inout), dimension(:, :), optional :: & + dqv_dt_prd, & ! Net production/loss effect + dqv_dt_qv ! water vapor effect integer :: nCol, nLev, iCol, iLev, iCf, kmax, kmin, k logical, dimension(size(p,1)) :: flg @@ -226,6 +229,15 @@ subroutine run(this, dt, p, h2opltc, h2o) h2o(iCol,iLev) = (h2oib(iCol) + (pltc(iCol,1)+pltc(iCol,3)*temp)*dt) / (1.0 + temp*dt) endif enddo + + if (present(dqv_dt_prd)) dqv_dt_prd(:, iLev) = pltc(:, 1) * dt + if (present(dqv_dt_qv)) then + where (pltc(:, 2) .ne. 0) + dqv_dt_qv(:, iLev) = (h2o(:, iLev) - pltc(:, 3)) * dt / pltc(:, 2) + elsewhere + dqv_dt_qv(:, iLev) = 0 + end where + end if enddo