From 34948e576921b556bd768e156c560742a4a81b45 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 23 Dec 2019 14:46:18 -0500 Subject: [PATCH 01/21] point sub-module to AnningCheng-NOAA/ccpp-physics MG3_v1 --- .gitmodules | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index 949d298df..9eb672592 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,4 +7,5 @@ url = https://github.com/NCAR/ccpp-framework [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/NCAR/ccpp-physics + url = https://github.com/AnningCheng-NOAA/ccpp-physics + branch = MG3_v1 From f2463184040b45e61444d4efaeb24b6ecf839f9f Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 23 Dec 2019 14:54:16 -0500 Subject: [PATCH 02/21] establish MG3_v1 branch in ufs weather model --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 904a43354..6478debf4 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 904a43354bc5922558dffc9b88acbc854c20e605 +Subproject commit 6478debf464519b50ea3bec29ca918f9f7ffba80 From 611e71b556ff340020c5d9f7e8be4119d6c738bf Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 23 Dec 2019 15:25:31 -0500 Subject: [PATCH 03/21] updated .gitmodules --- .gitmodules | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index 9eb672592..f8f14ba30 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,7 +4,8 @@ branch = dev/emc [submodule "ccpp/framework"] path = ccpp/framework - url = https://github.com/NCAR/ccpp-framework + url = https://github.com/AnningCheng-NOAA/ccpp-framework + branch = MG3_v1 [submodule "ccpp/physics"] path = ccpp/physics url = https://github.com/AnningCheng-NOAA/ccpp-physics From a240c2dd249381ab789c4b64ecb6803e31758e3d Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 31 Dec 2019 13:48:08 -0500 Subject: [PATCH 04/21] complete iccn changes in physics and GFS_layer --- atmos_cubed_sphere | 2 +- ccpp/physics | 2 +- gfsphysics/GFS_layer/GFS_driver.F90 | 12 +++---- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 2 +- gfsphysics/GFS_layer/GFS_typedefs.F90 | 38 ++++++++++----------- gfsphysics/GFS_layer/GFS_typedefs.meta | 10 +++--- gfsphysics/physics/m_micro_driver.F90 | 33 +++++++----------- gfsphysics/physics/micro_mg2_0.F90 | 11 +++--- gfsphysics/physics/micro_mg3_0.F90 | 11 +++--- 9 files changed, 57 insertions(+), 64 deletions(-) diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere index 452333aff..0e84f88b4 160000 --- a/atmos_cubed_sphere +++ b/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit 452333aff37f4c4348504bb9d485d98072e984e5 +Subproject commit 0e84f88b494b9e0a4097da50abe6b143330e8a2f diff --git a/ccpp/physics b/ccpp/physics index 6478debf4..829caff5c 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 6478debf464519b50ea3bec29ca918f9f7ffba80 +Subproject commit 829caff5c877b69c44fe94cdd6a33c86093143f7 diff --git a/gfsphysics/GFS_layer/GFS_driver.F90 b/gfsphysics/GFS_layer/GFS_driver.F90 index 3b6a94336..999df0f33 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -200,10 +200,10 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, & #ifndef CCPP call read_o3data (Model%ntoz, Model%me, Model%master) call read_h2odata (Model%h2o_phys, Model%me, Model%master) - if (Model%aero_in) then + if (Model%iaerclm) then call read_aerdata (Model%me,Model%master,Model%iflip,Model%idate) endif - if (Model%iccn) then + if (Model%iccn == 1) then call read_cidata ( Model%me, Model%master) endif #endif @@ -272,7 +272,7 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, & endif !--- read in and initialize IN and CCN - if (Model%iccn) then + if (Model%iccn == 1) then do nb = 1, nblks call setindxci (Init_parm%blksz(nb), Grid(nb)%xlat_d, Grid(nb)%jindx1_ci, & Grid(nb)%jindx2_ci, Grid(nb)%ddy_ci, Grid(nb)%xlon_d, & @@ -281,7 +281,7 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, & endif !--- read in and initialize aerosols - if (Model%aero_in) then + if (Model%iaerclm) then do nb = 1, nblks call setindxaer (Init_parm%blksz(nb),Grid(nb)%xlat_d,Grid(nb)%jindx1_aer, & Grid(nb)%jindx2_aer, Grid(nb)%ddy_aer, Grid(nb)%xlon_d, & @@ -959,7 +959,7 @@ subroutine GFS_phys_time_vary (Model, Grid, Tbd, Statein) endif !--- ICCN interpolation - if (Model%ICCN ) then + if (Model%ICCN == 1) then do nb = 1, nblks call ciinterpol (Model%me, blksz(nb), Model%idate, Model%fhour, & Grid(nb)%jindx1_ci, Grid(nb)%jindx2_ci, & @@ -971,7 +971,7 @@ subroutine GFS_phys_time_vary (Model, Grid, Tbd, Statein) endif !--- aerosol interpolation - if (Model%aero_in ) then + if (Model%iaerclm ) then do nb = 1, nblks call aerinterpol (Model%me, Model%master, blksz(nb), & Model%idate, Model%fhour, & diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 5b67f7faa..9bbd24a80 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -4691,7 +4691,7 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(1,1,2), Tbd%phy_f3d(1,1,3), & Tbd%phy_f3d(1,1,4), Tbd%phy_f3d(1,1,5), & Tbd%phy_f3d(1,1,kk), Tbd%aer_nm, & - Model%aero_in, Tbd%in_nm, Tbd%ccn_nm, Model%iccn, & + Tbd%in_nm, Tbd%ccn_nm, Model%iccn, & skip_macro, lprnt, & ! skip_macro, cn_prc, cn_snr, lprnt, & ! ipr, kdt, Grid%xlat, Grid%xlon) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 39520b0d4..73572f05c 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -595,7 +595,7 @@ module GFS_typedefs integer :: levrp1 !< number of vertical levels for radiation calculations plus one #endif integer :: nfxr !< second dimension for fluxr diagnostic variable (radiation) - logical :: aero_in !< flag for initializing aerosol data + logical :: iaerclm !< flag for initializing aerosol data #ifdef CCPP integer :: ntrcaer !< number of aerosol tracers for Morrison-Gettelman microphysics #endif @@ -1022,7 +1022,7 @@ module GFS_typedefs real(kind=kind_phys) :: julian !< current forecast julian date integer :: yearlen !< current length of the year ! - logical :: iccn !< using IN CCN forcing for MG2/3 + integer :: iccn !< using IN CCN forcing for MG2/3 #ifdef CCPP real(kind=kind_phys) :: sec !< seconds since model initialization real(kind=kind_phys), pointer :: si(:) !< vertical sigma coordinate for model initialization @@ -2653,8 +2653,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & real(kind=kind_phys) :: fhlwr = 3600. !< frequency for longwave radiation (secs) integer :: levr = -99 !< number of vertical levels for radiation calculations integer :: nfxr = 39+6 !< second dimension of input/output array fluxr - logical :: aero_in = .false. !< flag for initializing aero data - logical :: iccn = .false. !< logical to use IN CCN forcing for MG2/3 + logical :: iaerclm = .false. !< flag for initializing aero data + integer :: iccn = 0 !< logical to use IN CCN forcing for MG2/3 integer :: iflip = 1 !< iflip - is not the same as flipv integer :: isol = 0 !< use prescribed solar constant integer :: ico2 = 0 !< prescribed global mean value (old opernl) @@ -3012,7 +3012,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- coupling parameters cplflx, cplwav, cplchm, lsidea, & !--- radiation parameters - fhswr, fhlwr, levr, nfxr, aero_in, iflip, isol, ico2, ialb, & + fhswr, fhlwr, levr, nfxr, iaerclm, iflip, isol, ico2, ialb, & isot, iems, iaer, icliq_sw, iovr_sw, iovr_lw, ictm, isubc_sw,& isubc_lw, crick_proof, ccnorm, lwhtr, swhtr, & ! IN CCN forcing @@ -3224,17 +3224,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%levrp1 = Model%levr + 1 #endif Model%nfxr = nfxr - Model%aero_in = aero_in - if (Model%aero_in) then - ntrcaer = ntrcaerm - else - ntrcaer = 1 - endif -#ifdef CCPP - Model%ntrcaer = ntrcaer -#endif Model%iccn = iccn - if (Model%aero_in) Model%iccn = .false. ! further down: set Model%iccn to .false. ! for all microphysics schemes except ! MG2/3 (these are the only ones using ICCN) @@ -3244,6 +3234,15 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ialb = ialb Model%iems = iems Model%iaer = iaer + if (iaer/1000 == 1 .or. Model%iccn == 2) then + Model%iaerclm = .true. + ntrcaer = ntrcaerm + else + ntrcaer = 1 + endif +#ifdef CCPP + Model%ntrcaer = ntrcaer +#endif Model%icliq_sw = icliq_sw Model%iovr_sw = iovr_sw Model%iovr_lw = iovr_lw @@ -3269,7 +3268,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ncld = ncld Model%imp_physics = imp_physics ! turn off ICCN interpolation when MG2/3 are not used - if (.not. Model%imp_physics==Model%imp_physics_mg) Model%iccn = .false. + if (.not. Model%imp_physics==Model%imp_physics_mg) Model%iccn = 0 !--- Zhao-Carr MP parameters Model%psautco = psautco Model%prautco = prautco @@ -4074,7 +4073,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & endif if (Model%me == Model%master) & print *,' Using Morrison-Gettelman double moment microphysics', & - ' aero_in=', Model%aero_in, ' iccn=', Model%iccn, & + ' iaerclm=', Model%iaerclm, ' iccn=', Model%iccn, & ' mg_dcs=', Model%mg_dcs, ' mg_qcvar=', Model%mg_qcvar, & ' mg_ts_auto_ice=', Model%mg_ts_auto_ice, ' pdfflag=', Model%pdfflag, & ' mg_do_graupel=', Model%mg_do_graupel, ' mg_do_hail=', Model%mg_do_hail, & @@ -4276,7 +4275,6 @@ subroutine control_print(Model) print *, ' nslwr : ', Model%nslwr print *, ' levr : ', Model%levr print *, ' nfxr : ', Model%nfxr - print *, ' aero_in : ', Model%aero_in #ifdef CCPP print *, ' ntrcaer : ', Model%ntrcaer #endif @@ -4614,7 +4612,7 @@ subroutine grid_create (Grid, IM, Model) endif !--- iccn active - if ( Model%iccn ) then + if ( Model%iccn == 1) then allocate (Grid%ddy_ci (IM)) allocate (Grid%jindx1_ci (IM)) allocate (Grid%jindx2_ci (IM)) @@ -4624,7 +4622,7 @@ subroutine grid_create (Grid, IM, Model) endif !--- iaerclm active - if ( Model%aero_in ) then + if ( Model%iaerclm ) then allocate (Grid%ddy_aer (IM)) allocate (Grid%jindx1_aer(IM)) allocate (Grid%jindx2_aer(IM)) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.meta b/gfsphysics/GFS_layer/GFS_typedefs.meta index 79ab00129..b62120bde 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -2030,9 +2030,9 @@ units = count dimensions = () type = integer -[aero_in] - standard_name = flag_for_aerosol_input_MG - long_name = flag for using aerosols in Morrison-Gettelman MP +[iaerclm] + standard_name = flag_for_aerosol_input_MG_radiation + long_name = flag for using aerosols in Morrison-Gettelman MP_radiation units = flag dimensions = () type = logical @@ -3636,9 +3636,9 @@ [iccn] standard_name = flag_for_in_ccn_forcing_for_morrison_gettelman_microphysics long_name = flag for IN and CCN forcing for morrison gettelman microphysics - units = flag + units = none dimensions = () - type = logical + type = integer [sec] standard_name = seconds_elapsed_since_model_initialization long_name = seconds elapsed since model initialization diff --git a/gfsphysics/physics/m_micro_driver.F90 b/gfsphysics/physics/m_micro_driver.F90 index 9d4d3a318..715c8c9bc 100644 --- a/gfsphysics/physics/m_micro_driver.F90 +++ b/gfsphysics/physics/m_micro_driver.F90 @@ -14,7 +14,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & &, CLLS_io, KCBL & &, CLDREFFL, CLDREFFI, CLDREFFR, CLDREFFS & &, CLDREFFG, aerfld_i & - &, aero_in, naai_i, npccn_i, iccn & + &, naai_i, npccn_i, iccn & &, skip_macro & ! &, skip_macro, cn_prc2, cn_snr & &, lprnt, alf_fac, qc_min, pdfflag & @@ -62,7 +62,9 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & integer, parameter :: ncolmicro = 1 integer,intent(in) :: im, ix,lm, ipr, kdt, fprcp, pdfflag - logical,intent(in) :: flipv, aero_in, skip_macro, lprnt, iccn + logical,intent(in) :: flipv, skip_macro, lprnt + integer,intent(in) : + integer,intent(in) ::: iccn real (kind=kind_phys), intent(in):: dt_i, alf_fac, qc_min(2) real (kind=kind_phys), dimension(ix,lm),intent(in) :: & @@ -304,7 +306,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,ll) radheat(i,k) = lwheat_i(i,ll) + swheat_i(i,ll) rhc(i,k) = rhc_i(i,ll) - if (iccn) then + if (iccn == 1) then CDNC_NUC(i,k) = npccn_i(i,ll) INC_NUC(i,k) = naai_i (i,ll) endif @@ -365,7 +367,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,k) radheat(i,k) = lwheat_i(i,k) + swheat_i(i,k) rhc(i,k) = rhc_i(i,k) - if (iccn) then + if (iccn == 1) then CDNC_NUC(i,k) = npccn_i(i,k) INC_NUC(i,k) = naai_i (i,k) endif @@ -528,7 +530,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & NCPL(i,l) = MAX( NCPL(i,l), 0.) NCPI(i,l) = MAX( NCPI(i,l), 0.) RAD_CF(i,l) = max(0.0, min(CLLS(i,l)+CLCN(i,l), 1.0)) - if (.not. iccn) then + if (iccn.ne.1) then CDNC_NUC(i,l) = 0.0 INC_NUC(i,l) = 0.0 endif @@ -583,7 +585,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! allocate(AERMASSMIX(IM,LM,15)) - if ( aero_in ) then + if (iccn == 2) then AERMASSMIX(:,:,1:ntrcaer) = aerfld_i(:,:,1:ntrcaer) else AERMASSMIX(:,:,1:5) = 1.e-6 @@ -766,12 +768,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & tauxr8 = ter8(K) endif -! if(aero_in) then AeroAux = AeroProps(I, K) -! else -! call init_Aer(AeroAux) -! call init_Aer(AeroAux_b) -! endif pfrz_inc_r8(k) = 0.0 rh1_r8 = 0.0 !related to cnv_dql_dt, needed to changed soon @@ -846,7 +843,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & SC_ICE(i,k) = ((tice-temp(i,k))*tx1 + (temp(i,k)-t_ice_all)*rhc(i,k)) & * t_ice_denom endif - if (.not. iccn) then + if (iccn.ne.1) then CDNC_NUC(I,k) = npccninr8(k) INC_NUC (I,k) = naair8(k) endif @@ -981,7 +978,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! temp(i,k) = th1(i,k) * PK(i,k) RAD_CF(i,k) = min(CLLS(i,k)+CLCN(i,k), 1.0) ! - if (.not. iccn) then + if (iccn.ne.1) then if (PFRZ(i,k) > 0.0) then INC_NUC(i,k) = INC_NUC(i,k) * PFRZ(i,k) NHET_NUC(i,k) = NHET_NUC(i,k) * PFRZ(i,k) @@ -1130,11 +1127,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & endif -! if(aero_in) then - AeroAux = AeroProps(I, K) -! else -! call init_Aer(AeroAux) -! end if + AeroAux = AeroProps(I, K) call getINsubset(1, AeroAux, AeroAux_b) naux = AeroAux_b%nmods if (nbincontactdust < naux) then @@ -1348,7 +1341,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & drout2, dsout2, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, & & lev_sed_strt) ! LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) @@ -1482,7 +1475,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & qgout2, ngout2, dgout2, freqg, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, & & lev_sed_strt) LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) diff --git a/gfsphysics/physics/micro_mg2_0.F90 b/gfsphysics/physics/micro_mg2_0.F90 index 325a2dbbe..82d34ed9f 100644 --- a/gfsphysics/physics/micro_mg2_0.F90 +++ b/gfsphysics/physics/micro_mg2_0.F90 @@ -403,7 +403,7 @@ subroutine micro_mg_tend ( & drout2, dsout2, & freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) + prer_evap, xlat, xlon, lprnt, iccn, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -474,7 +474,8 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt, iccn, aero_in + logical, intent(in) :: lprnt + integer, intent(in) :: iccn ! used for scavenging @@ -1123,7 +1124,7 @@ subroutine micro_mg_tend ( & enddo enddo - if(iccn) then + if(iccn == 1) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1162,7 +1163,7 @@ subroutine micro_mg_tend ( & ncal = zero end where - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1177,7 +1178,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (aero_in) then + elseif (iccn == 2) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then diff --git a/gfsphysics/physics/micro_mg3_0.F90 b/gfsphysics/physics/micro_mg3_0.F90 index cbd25370a..db0870bee 100644 --- a/gfsphysics/physics/micro_mg3_0.F90 +++ b/gfsphysics/physics/micro_mg3_0.F90 @@ -505,7 +505,7 @@ subroutine micro_mg_tend ( & !--ag freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) + prer_evap, xlat, xlon, lprnt, iccn, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -593,7 +593,8 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt, iccn, aero_in + logical, intent(in) :: lprnt + integer, intent(in) :: iccn ! used for scavenging @@ -1441,7 +1442,7 @@ subroutine micro_mg_tend ( & enddo enddo ! - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1495,7 +1496,7 @@ subroutine micro_mg_tend ( & enddo enddo - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1510,7 +1511,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (aero_in) then + elseif (iccn == 2) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then From ea47b607007d3621d3da42d80bee492bbe23af3b Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Thu, 2 Jan 2020 14:06:11 -0500 Subject: [PATCH 05/21] passed ccpp compilation and testing ipd compilation --- .gitmodules | 6 +- atmos_model.F90 | 34 +++++- ccpp/physics | 2 +- fv3_cap.F90 | 18 ++- gfsphysics/GFS_layer/GFS_driver.F90 | 2 +- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 8 +- gfsphysics/GFS_layer/GFS_radiation_driver.F90 | 38 ++++-- gfsphysics/physics/m_micro_driver.F90 | 33 +++--- gfsphysics/physics/micro_mg2_0.F90 | 11 +- gfsphysics/physics/micro_mg3_0.F90 | 11 +- gfsphysics/physics/sflx.f | 1 + gfsphysics/physics/ugwp_driver_v0.f | 4 +- io/FV3GFS_io.F90 | 4 +- io/post_gfs.F90 | 81 ++++++++++--- module_fcst_grid_comp.F90 | 109 ++++++++++++------ module_fv3_config.F90 | 7 +- 16 files changed, 261 insertions(+), 108 deletions(-) diff --git a/.gitmodules b/.gitmodules index f8f14ba30..949d298df 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,9 +4,7 @@ branch = dev/emc [submodule "ccpp/framework"] path = ccpp/framework - url = https://github.com/AnningCheng-NOAA/ccpp-framework - branch = MG3_v1 + url = https://github.com/NCAR/ccpp-framework [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/AnningCheng-NOAA/ccpp-physics - branch = MG3_v1 + url = https://github.com/NCAR/ccpp-physics diff --git a/atmos_model.F90 b/atmos_model.F90 index 23e30e76c..34e1a6c99 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -113,7 +113,8 @@ module atmos_model_mod FV3GFS_diag_register, FV3GFS_diag_output, & DIAG_SIZE use fv_iau_mod, only: iau_external_data_type,getiauforcing,iau_initialize -use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout +use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout, & + frestart, restart_endfcst !----------------------------------------------------------------------- @@ -221,7 +222,8 @@ module atmos_model_mod logical,parameter :: flip_vc = .true. #endif - real(kind=IPD_kind_phys), parameter :: zero=0.0, one=1.0 + real(kind=IPD_kind_phys), parameter :: zero = 0.0_IPD_kind_phys, & + one = 1.0_IPD_kind_phys contains @@ -944,7 +946,7 @@ end subroutine update_atmos_model_state subroutine atmos_model_end (Atmos) type (atmos_data_type), intent(inout) :: Atmos !---local variables - integer :: idx + integer :: idx, seconds #ifdef CCPP integer :: ierr #endif @@ -952,9 +954,11 @@ subroutine atmos_model_end (Atmos) !----------------------------------------------------------------------- !---- termination routine for atmospheric model ---- - call atmosphere_end (Atmos % Time, Atmos%grid) - call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, & - IPD_Control, Atmos%domain) + call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst) + if(restart_endfcst) then + call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, & + IPD_Control, Atmos%domain) + endif #ifdef CCPP ! Fast physics (from dynamics) are finalized in atmosphere_end above; @@ -1457,6 +1461,24 @@ subroutine update_atmos_chemistry(state, rc) enddo enddo + ! -- zero out accumulated fields +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, IPD_Control, IPD_Data) & +!$OMP private (j, jb, i, ib, nb, ix) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%coupling%rainc_cpl(ix) = zero + if (.not.IPD_Control%cplflx) then + IPD_Data(nb)%coupling%rain_cpl(ix) = zero + IPD_Data(nb)%coupling%snow_cpl(ix) = zero + end if + enddo + enddo + if (IPD_Control%debug) then ! -- diagnostics write(6,'("update_atmos: prsi - min/max/avg",3g16.6)') minval(prsi), maxval(prsi), sum(prsi)/size(prsi) diff --git a/ccpp/physics b/ccpp/physics index 829caff5c..b30d0ce83 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 829caff5c877b69c44fe94cdd6a33c86093143f7 +Subproject commit b30d0ce83f25e80616770dd433451b99eccb3a57 diff --git a/fv3_cap.F90 b/fv3_cap.F90 index d04abdce2..a099ae32e 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -30,7 +30,7 @@ module fv3gfs_cap_mod calendar, calendar_type, cpl, & force_date_from_configure, & cplprint_flag,output_1st_tstep_rst, & - first_kdt + first_kdt,num_restart_interval use module_fv3_io_def, only: num_pes_fcst,write_groups,app_domain, & num_files, filename_base, & @@ -278,9 +278,16 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) CALL ESMF_ConfigLoadFile(config=CF ,filename='model_configure' ,rc=RC) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! - CALL ESMF_ConfigGetAttribute(config=CF,value=restart_interval, & - label ='restart_interval:',rc=rc) + num_restart_interval = ESMF_ConfigGetLen(config=CF, label ='restart_interval:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if(mype == 0) print *,'af nems config,num_restart_interval=',num_restart_interval + if (num_restart_interval<=0) num_restart_interval = 1 + allocate(restart_interval(num_restart_interval)) + restart_interval = 0 + CALL ESMF_ConfigGetAttribute(CF,valueList=restart_interval,label='restart_interval:', & + count=num_restart_interval, rc=RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if(mype == 0) print *,'af nems config,restart_interval=',restart_interval ! CALL ESMF_ConfigGetAttribute(config=CF,value=calendar, & label ='calendar:',rc=rc) @@ -326,9 +333,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) label ='app_domain:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - if(mype == 0) print *,'af nems config,restart_interval=',restart_interval, & - 'quilting=',quilting,'write_groups=',write_groups,wrttasks_per_group, & - 'calendar=',trim(calendar),'calendar_type=',calendar_type + if(mype == 0) print *,'af nems config,quilting=',quilting,'write_groups=', & + write_groups,wrttasks_per_group,'calendar=',trim(calendar),'calendar_type=',calendar_type ! CALL ESMF_ConfigGetAttribute(config=CF,value=num_files, & label ='num_files:',rc=rc) diff --git a/gfsphysics/GFS_layer/GFS_driver.F90 b/gfsphysics/GFS_layer/GFS_driver.F90 index 999df0f33..644816ea9 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -422,7 +422,7 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, & call cires_ugwp_init(Model%me, Model%master, Model%nlunit, Init_parm%logunit, & Model%fn_nml, Model%lonr, Model%latr, Model%levs, & Init_parm%ak, Init_parm%bk, p_ref, Model%dtp, & - Model%cdmbgwd, Model%cgwf, Model%prslrd0, Model%ral_ts) + Model%cdmbgwd(1:2), Model%cgwf, Model%prslrd0, Model%ral_ts) endif #endif diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 9bbd24a80..b728a7237 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -3184,9 +3184,13 @@ subroutine GFS_physics_driver & dtdt(1:im,:) = Stateout%gt0(1:im,:) endif ! end if_ldiag3d/cnvgwd - if (Model%ldiag3d) then + if (Model%ldiag3d .or. Model%cplchm) then dqdt(1:im,:,1) = Stateout%gq0(1:im,:,1) - endif ! end if_ldiag3d + endif ! end if_ldiag3d/cplchm + + if (Model%cplchm) then + Coupling%dqdti(1:im,:) = zero + endif ! end if_cplchm #ifdef GFS_HYDRO call get_phi(im, ix, levs, ntrac, Stateout%gt0, Stateout%gq0, & diff --git a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 index c7323d6bb..5215f1b2b 100644 --- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 @@ -2116,18 +2116,40 @@ subroutine GFS_radiation_driver & Diag%fluxr(i,11-j) = Diag%fluxr(i,11-j) + tem0d * Statein%prsi(i,itop+kt) Diag%fluxr(i,14-j) = Diag%fluxr(i,14-j) + tem0d * Statein%prsi(i,ibtc+kb) Diag%fluxr(i,17-j) = Diag%fluxr(i,17-j) + tem0d * Statein%tgrs(i,itop) + enddo + enddo ! Anning adds optical depth and emissivity output - tem1 = 0. - tem2 = 0. - do k=ibtc,itop - tem1 = tem1 + cldtausw(i,k) ! approx .55 mu channel - tem2 = tem2 + cldtaulw(i,k) ! approx 10. mu channel + if (Model%lsswr .and. (nday > 0)) then + do j = 1, 3 + do i = 1, IM + tem0d = raddt * cldsa(i,j) + itop = mtopa(i,j) - kd + ibtc = mbota(i,j) - kd + tem1 = 0. + do k=ibtc,itop + tem1 = tem1 + cldtausw(i,k) ! approx .55 um channel + enddo + Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1 enddo - Diag%fluxr(i,43-j) = Diag%fluxr(i,43-j) + tem0d * tem1 - Diag%fluxr(i,46-j) = Diag%fluxr(i,46-j) + tem0d * (1.0-exp(-tem2)) enddo - enddo + endif + + if (Model%lslwr) then + do j = 1, 3 + do i = 1, IM + tem0d = raddt * cldsa(i,j) + itop = mtopa(i,j) - kd + ibtc = mbota(i,j) - kd + tem2 = 0. + do k=ibtc,itop + tem2 = tem2 + cldtaulw(i,k) ! approx 10. um channel + enddo + Diag%fluxr(i,46-j) = Diag%fluxr(i,46-j) + tem0d * (1.0-exp(-tem2)) + enddo + enddo + endif + endif endif ! end_if_lssav diff --git a/gfsphysics/physics/m_micro_driver.F90 b/gfsphysics/physics/m_micro_driver.F90 index 715c8c9bc..9d4d3a318 100644 --- a/gfsphysics/physics/m_micro_driver.F90 +++ b/gfsphysics/physics/m_micro_driver.F90 @@ -14,7 +14,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & &, CLLS_io, KCBL & &, CLDREFFL, CLDREFFI, CLDREFFR, CLDREFFS & &, CLDREFFG, aerfld_i & - &, naai_i, npccn_i, iccn & + &, aero_in, naai_i, npccn_i, iccn & &, skip_macro & ! &, skip_macro, cn_prc2, cn_snr & &, lprnt, alf_fac, qc_min, pdfflag & @@ -62,9 +62,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & integer, parameter :: ncolmicro = 1 integer,intent(in) :: im, ix,lm, ipr, kdt, fprcp, pdfflag - logical,intent(in) :: flipv, skip_macro, lprnt - integer,intent(in) : - integer,intent(in) ::: iccn + logical,intent(in) :: flipv, aero_in, skip_macro, lprnt, iccn real (kind=kind_phys), intent(in):: dt_i, alf_fac, qc_min(2) real (kind=kind_phys), dimension(ix,lm),intent(in) :: & @@ -306,7 +304,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,ll) radheat(i,k) = lwheat_i(i,ll) + swheat_i(i,ll) rhc(i,k) = rhc_i(i,ll) - if (iccn == 1) then + if (iccn) then CDNC_NUC(i,k) = npccn_i(i,ll) INC_NUC(i,k) = naai_i (i,ll) endif @@ -367,7 +365,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,k) radheat(i,k) = lwheat_i(i,k) + swheat_i(i,k) rhc(i,k) = rhc_i(i,k) - if (iccn == 1) then + if (iccn) then CDNC_NUC(i,k) = npccn_i(i,k) INC_NUC(i,k) = naai_i (i,k) endif @@ -530,7 +528,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & NCPL(i,l) = MAX( NCPL(i,l), 0.) NCPI(i,l) = MAX( NCPI(i,l), 0.) RAD_CF(i,l) = max(0.0, min(CLLS(i,l)+CLCN(i,l), 1.0)) - if (iccn.ne.1) then + if (.not. iccn) then CDNC_NUC(i,l) = 0.0 INC_NUC(i,l) = 0.0 endif @@ -585,7 +583,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! allocate(AERMASSMIX(IM,LM,15)) - if (iccn == 2) then + if ( aero_in ) then AERMASSMIX(:,:,1:ntrcaer) = aerfld_i(:,:,1:ntrcaer) else AERMASSMIX(:,:,1:5) = 1.e-6 @@ -768,7 +766,12 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & tauxr8 = ter8(K) endif +! if(aero_in) then AeroAux = AeroProps(I, K) +! else +! call init_Aer(AeroAux) +! call init_Aer(AeroAux_b) +! endif pfrz_inc_r8(k) = 0.0 rh1_r8 = 0.0 !related to cnv_dql_dt, needed to changed soon @@ -843,7 +846,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & SC_ICE(i,k) = ((tice-temp(i,k))*tx1 + (temp(i,k)-t_ice_all)*rhc(i,k)) & * t_ice_denom endif - if (iccn.ne.1) then + if (.not. iccn) then CDNC_NUC(I,k) = npccninr8(k) INC_NUC (I,k) = naair8(k) endif @@ -978,7 +981,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! temp(i,k) = th1(i,k) * PK(i,k) RAD_CF(i,k) = min(CLLS(i,k)+CLCN(i,k), 1.0) ! - if (iccn.ne.1) then + if (.not. iccn) then if (PFRZ(i,k) > 0.0) then INC_NUC(i,k) = INC_NUC(i,k) * PFRZ(i,k) NHET_NUC(i,k) = NHET_NUC(i,k) * PFRZ(i,k) @@ -1127,7 +1130,11 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & endif - AeroAux = AeroProps(I, K) +! if(aero_in) then + AeroAux = AeroProps(I, K) +! else +! call init_Aer(AeroAux) +! end if call getINsubset(1, AeroAux, AeroAux_b) naux = AeroAux_b%nmods if (nbincontactdust < naux) then @@ -1341,7 +1348,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & drout2, dsout2, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & & lev_sed_strt) ! LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) @@ -1475,7 +1482,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & qgout2, ngout2, dgout2, freqg, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & & lev_sed_strt) LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) diff --git a/gfsphysics/physics/micro_mg2_0.F90 b/gfsphysics/physics/micro_mg2_0.F90 index 82d34ed9f..325a2dbbe 100644 --- a/gfsphysics/physics/micro_mg2_0.F90 +++ b/gfsphysics/physics/micro_mg2_0.F90 @@ -403,7 +403,7 @@ subroutine micro_mg_tend ( & drout2, dsout2, & freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, nlball) + prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -474,8 +474,7 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt - integer, intent(in) :: iccn + logical, intent(in) :: lprnt, iccn, aero_in ! used for scavenging @@ -1124,7 +1123,7 @@ subroutine micro_mg_tend ( & enddo enddo - if(iccn == 1) then + if(iccn) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1163,7 +1162,7 @@ subroutine micro_mg_tend ( & ncal = zero end where - if (iccn == 1) then + if (iccn) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1178,7 +1177,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (iccn == 2) then + elseif (aero_in) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then diff --git a/gfsphysics/physics/micro_mg3_0.F90 b/gfsphysics/physics/micro_mg3_0.F90 index db0870bee..cbd25370a 100644 --- a/gfsphysics/physics/micro_mg3_0.F90 +++ b/gfsphysics/physics/micro_mg3_0.F90 @@ -505,7 +505,7 @@ subroutine micro_mg_tend ( & !--ag freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, nlball) + prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -593,8 +593,7 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt - integer, intent(in) :: iccn + logical, intent(in) :: lprnt, iccn, aero_in ! used for scavenging @@ -1442,7 +1441,7 @@ subroutine micro_mg_tend ( & enddo enddo ! - if (iccn == 1) then + if (iccn) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1496,7 +1495,7 @@ subroutine micro_mg_tend ( & enddo enddo - if (iccn == 1) then + if (iccn) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1511,7 +1510,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (iccn == 2) then + elseif (aero_in) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then diff --git a/gfsphysics/physics/sflx.f b/gfsphysics/physics/sflx.f index 29467fe37..f84006be9 100644 --- a/gfsphysics/physics/sflx.f +++ b/gfsphysics/physics/sflx.f @@ -251,6 +251,7 @@ subroutine sflx & runoff2 = 0.0 runoff3 = 0.0 snomlt = 0.0 + rc = 0.0 ! --- ... define local variable ice to achieve: ! sea-ice case, ice = 1 diff --git a/gfsphysics/physics/ugwp_driver_v0.f b/gfsphysics/physics/ugwp_driver_v0.f index 804bbac19..0ca37e818 100644 --- a/gfsphysics/physics/ugwp_driver_v0.f +++ b/gfsphysics/physics/ugwp_driver_v0.f @@ -46,7 +46,9 @@ subroutine cires_ugwp_driver_v0(me, master, &, rain real(kind=kind_phys), intent(in), dimension(im,levs) :: ugrs - &, vgrs, tgrs, qgrs, prsi, prsl, prslk, phii, phil, del + &, vgrs, tgrs, qgrs, prsl, prslk, phil, del + real(kind=kind_phys), intent(in), dimension(im,levs+1) :: prsi + &, phii ! real(kind=kind_phys), intent(in) :: oro_stat(im,nmtvr) real(kind=kind_phys), intent(in), dimension(im) :: hprime, oc diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 4b2d1e426..e8c980f4d 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -157,10 +157,10 @@ subroutine FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, Model, fv_dom type(domain2d), intent(in) :: fv_domain character(len=32), optional, intent(in) :: timestamp - !--- read in surface data from chgres + !--- write surface data from chgres call sfc_prop_restart_write (IPD_Data%Sfcprop, Atm_block, Model, fv_domain, timestamp) - !--- read in physics restart data + !--- write physics restart data call phys_restart_write (IPD_Restart, Atm_block, Model, fv_domain, timestamp) end subroutine FV3GFS_restart_write diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index 9f3a0a80f..98d4fef50 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -12,6 +12,7 @@ module post_gfs include 'mpif.h' integer mype, nbdl + logical setvar_atmfile, setvar_sfcfile, read_postcntrl public post_run_gfs, post_getattr_gfs contains @@ -28,9 +29,10 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & ! use ctlblk_mod, only : komax,ifhr,ifmin,modelname,datapd,fld_info, & npset,grib,gocart_on,icount_calmict, jsta, & - jend,im, nsoil + jend,im, nsoil, filenameflat use gridspec_mod, only : maptype, gridtype use grib2_module, only : gribit2,num_pset,nrecout,first_grbtbl + use xml_perl_data,only : paramset ! !----------------------------------------------------------------------- ! @@ -53,9 +55,8 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & integer n,nwtpg,ieof,lcntrl,ierr,i,j,k,jts,jte,mynsoil integer,allocatable :: jstagrp(:),jendgrp(:) integer,save :: kpo,kth,kpv + logical,save :: log_postalct=.false. real,dimension(komax),save :: po, th, pv - logical,save :: log_postalct=.false. - logical,save :: setvar_atmfile=.false.,setvar_sfcfile=.false. logical :: Log_runpost character(255) :: post_fname*255 @@ -124,6 +125,7 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & ! log_postalct = .true. first_grbtbl = .true. + read_postcntrl = .true. ! ENDIF ! @@ -135,6 +137,8 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & ifmin = mynfmin if (ifhr == 0 ) ifmin = 0 if(mype==0) print *,'bf set_postvars,ifmin=',ifmin,'ifhr=',ifhr + setvar_atmfile=.false. + setvar_sfcfile=.false. call set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & setvar_sfcfile) @@ -145,8 +149,28 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & ! 20190807 no need to call microinit for GFDLMP ! call MICROINIT ! - if(grib=="grib2" .and. first_grbtbl) then - call read_xml() + if(grib=="grib2" .and. read_postcntrl) then + if (ifhr == 0) then + filenameflat = 'postxconfig-NT_FH00.txt' + call read_xml() + if(mype==0) print *,'af read_xml at fh00,name=',trim(filenameflat) + else if(ifhr > 0) then + filenameflat = 'postxconfig-NT.txt' + if(size(paramset)>0) then + do i=1,size(paramset) + if (size(paramset(i)%param)>0) then + deallocate(paramset(i)%param) + nullify(paramset(i)%param) + endif + enddo + deallocate(paramset) + nullify(paramset) + endif + num_pset = 0 + call read_xml() + if(mype==0) print *,'af read_xml,name=',trim(filenameflat),'ifhr=',ifhr + read_postcntrl = .false. + endif endif ! IEOF = 0 @@ -181,9 +205,6 @@ subroutine post_run_gfs(wrt_int_state,mypei,mpicomp,lead_write, & endif ! enddo -! - setvar_atmfile = .false. - setvar_sfcfile = .false. ! endif @@ -335,7 +356,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & avgetrans, avgesnow, avgprec_cont, avgcprate_cont,& avisbeamswin, avisdiffswin, airbeamswin, airdiffswin, & alwoutc, alwtoac, aswoutc, aswtoac, alwinc, aswinc,& - avgpotevp, snoavg, si, cuppt + avgpotevp, snoavg, ti, si, cuppt use soil, only: sldpth, sh2o, smc, stc use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice use ctlblk_mod, only: im, jm, lm, lp1, jsta, jend, jsta_2l, jend_2u, jsta_m,jend_m, & @@ -477,6 +498,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & qs(i,j) = SPVAL twbs(i,j) = SPVAL qwbs(i,j) = SPVAL + ths(i,j) = SPVAL enddo enddo @@ -1122,9 +1144,30 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & !$omp parallel do private(i,j) do j=jsta,jend do i=ista, iend - sr(i,j) = arrayr42d(i,j) + if (arrayr42d(i,j) /= spval) then + !set range within (0,1) + sr(i,j) = min(1.,max(0.,sr(i,j))) + else + sr(i,j) = spval + endif + enddo + enddo + endif + + ! sea ice skin temperature + if(trim(fieldname)=='tisfc') then + !$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (arrayr42d(i,j) /= spval) then + ti(i,j) = arrayr42d(i,j) + if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval + else + ti(i,j) = spval + endif enddo enddo +! print *,'in gfs_post, get tisfc=',maxval(ti), minval(ti) endif ! vegetation fraction @@ -1237,7 +1280,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do j=jsta,jend do i=ista, iend stc(i,j,1) = arrayr42d(i,j) - if (sm(i,j) /= 0.0) stc(i,j,1) = spval + !mask open water areas, combine with sea ice tmp + if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval enddo enddo endif @@ -1248,7 +1292,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do j=jsta,jend do i=ista, iend stc(i,j,2) = arrayr42d(i,j) - if (sm(i,j) /= 0.0) stc(i,j,2) = spval + !mask open water areas, combine with sea ice tmp + if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval enddo enddo endif @@ -1259,7 +1304,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do j=jsta,jend do i=ista, iend stc(i,j,3) = arrayr42d(i,j) - if (sm(i,j) /= 0.0) stc(i,j,3) = spval + !mask open water areas, combine with sea ice tmp + if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval enddo enddo endif @@ -1270,7 +1316,8 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do j=jsta,jend do i=ista, iend stc(i,j,4) = arrayr42d(i,j) - if (sm(i,j) /= 0.0) stc(i,j,4) = spval + !mask open water areas, combine with sea ice tmp + if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval enddo enddo endif @@ -2313,6 +2360,12 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & !$omp parallel do private(i,j) do j=jsta,jend do i=ista, iend + !assign sst + if (sm(i,j) /= 0.0 .and. ths(i,j) /= spval) then + sst(i,j) = ths(i,j) + else + sst(i,j) = spval + endif if (ths(i,j) /= spval) then ths(i,j) = ths(i,j)* (p1000/pint(i,j,lp1))**capa thz0(i,j) = ths(i,j) diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90 index 6ff1f39c7..fef9698ab 100644 --- a/module_fcst_grid_comp.F90 +++ b/module_fcst_grid_comp.F90 @@ -22,11 +22,12 @@ module module_fcst_grid_comp ! use time_manager_mod, only: time_type, set_calendar_type, set_time, & set_date, days_in_month, month_name, & - operator(+), operator (<), operator (>), & - operator (/=), operator (/), operator (==),& - operator (*), THIRTY_DAY_MONTHS, JULIAN, & - NOLEAP, NO_CALENDAR, date_to_string, & - get_date + operator(+), operator(-), operator (<), & + operator (>), operator (/=), operator (/), & + operator (==), operator (*), & + THIRTY_DAY_MONTHS, JULIAN, NOLEAP, & + NO_CALENDAR, date_to_string, get_date, & + get_time use atmos_model_mod, only: atmos_model_init, atmos_model_end, & get_atmos_model_ungridded_dim, & @@ -70,7 +71,8 @@ module module_fcst_grid_comp iau_offset use module_fv3_config, only: dt_atmos, calendar, restart_interval, & quilting, calendar_type, cpl, & - cplprint_flag, force_date_from_configure + cplprint_flag, force_date_from_configure, & + num_restart_interval, frestart, restart_endfcst ! !----------------------------------------------------------------------- ! @@ -88,7 +90,8 @@ module module_fcst_grid_comp type(atmos_data_type) :: Atm type(time_type) :: Time_atmos, Time_init, Time_end, & Time_step_atmos, Time_step_ocean, & - Time_restart, Time_step_restart + Time_restart, Time_step_restart, & + Time_atstart integer :: num_atmos_calls, ret, intrm_rst end type @@ -179,12 +182,11 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) integer :: Run_length integer,dimension(6) :: date, date_end - integer :: res_intvl integer :: mpi_comm_comp ! logical,save :: first=.true. character(len=9) :: month - integer :: initClock, unit, nfhour + integer :: initClock, unit, nfhour, total_inttime integer :: mype, ntasks character(3) cfhour character(4) dateSY @@ -203,7 +205,8 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) real(ESMF_KIND_R8),parameter :: dtor = 180.0_ESMF_KIND_R8 / 3.1415926535897931_ESMF_KIND_R8 integer :: jsc, jec, isc, iec, nlev type(domain2D) :: domain - integer :: n, fcstNpes + integer :: n, fcstNpes, tmpvar + logical :: single_restart integer, allocatable, dimension(:) :: isl, iel, jsl, jel integer, allocatable, dimension(:,:,:) :: deBlockList @@ -271,8 +274,9 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) !----------------------------------------------------------------------- ! call ESMF_ClockGet(clock, CurrTime=CurrTime, StartTime=StartTime, & - StopTime=StopTime, RunDuration=RunDuration, rc=rc) + StopTime=StopTime, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + RunDuration = StopTime - CurrTime date_init = 0 call ESMF_TimeGet (StartTime, & @@ -317,16 +321,46 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) ! atm_int_state%Time_step_atmos = set_time (dt_atmos,0) atm_int_state%num_atmos_calls = Run_length / dt_atmos + atm_int_state%Time_atstart = atm_int_state%Time_atmos if (mype == 0) write(0,*)'num_atmos_calls=',atm_int_state%num_atmos_calls,'time_init=', & date_init,'time_atmos=',date,'time_end=',date_end,'dt_atmos=',dt_atmos, & 'Run_length=',Run_length - res_intvl = restart_interval*3600 - atm_int_state%Time_step_restart = set_time (res_intvl, 0) - atm_int_state%Time_restart = atm_int_state%Time_atmos + atm_int_state%Time_step_restart - atm_int_state%intrm_rst = 0 - if (res_intvl>0) atm_int_state%intrm_rst = 1 - atm_int_state%Atm%iau_offset = iau_offset -! + frestart = 0 + single_restart = .false. + call get_time(atm_int_state%Time_end - atm_int_state%Time_atstart,total_inttime) + if(num_restart_interval == 2) then + if(restart_interval(2)== -1) single_restart = .true. + endif + if(single_restart) then + frestart(1) = restart_interval(1) * 3600 + elseif ( num_restart_interval == 1) then + if(restart_interval(1) == 0) then + frestart(1) = total_inttime + else if(restart_interval(1) > 0) then + tmpvar = restart_interval(1) * 3600 + frestart(1) = tmpvar + atm_int_state%Time_step_restart = set_time (tmpvar, 0) + atm_int_state%Time_restart = atm_int_state%Time_atstart + atm_int_state%Time_step_restart + i = 2 + do while ( atm_int_state%Time_restart < atm_int_state%Time_end ) + frestart(i) = frestart(i-1) + tmpvar + atm_int_state%Time_restart = atm_int_state%Time_restart + atm_int_state%Time_step_restart + i = i + 1 + enddo + endif + else if(num_restart_interval > 1) then + do i=1,num_restart_interval + frestart(i) = restart_interval(i) * 3600 + enddo + endif + restart_endfcst = .false. + if ( ANY(frestart(:) == total_inttime) ) restart_endfcst = .true. + if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'restart_endfcst=',restart_endfcst, & + 'total_inttime=',total_inttime + + atm_int_state%intrm_rst = 0 + if (frestart(1)>0) atm_int_state%intrm_rst = 1 + atm_int_state%Atm%iau_offset = iau_offset ! !----- write time stamps (for start time and end time) ------ @@ -737,9 +771,10 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) !----------------------------------------------------------------------- !*** local variables ! - integer :: i,j, mype, na, date(6) + integer :: i,j, mype, na, date(6), seconds character(20) :: compname - + + type(time_type) :: restart_inctime type(ESMF_Time) :: currtime integer(kind=ESMF_KIND_I8) :: ntimestep_esmf character(len=64) :: timestamp @@ -776,13 +811,16 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) !--- intermediate restart if (atm_int_state%intrm_rst>0) then - if ((na /= atm_int_state%num_atmos_calls) .and. & - (atm_int_state%Time_atmos == atm_int_state%Time_restart)) then - timestamp = date_to_string (atm_int_state%Time_restart) - call atmos_model_restart(atm_int_state%Atm, timestamp) - - call wrt_atmres_timestamp(atm_int_state,timestamp) - atm_int_state%Time_restart = atm_int_state%Time_restart + atm_int_state%Time_step_restart + if (na /= atm_int_state%num_atmos_calls-1) then + call get_time(atm_int_state%Time_atmos - atm_int_state%Time_atstart, seconds) + if (ANY(frestart(:) == seconds)) then + restart_inctime = set_time(seconds, 0) + atm_int_state%Time_restart = atm_int_state%Time_atstart + restart_inctime + timestamp = date_to_string (atm_int_state%Time_restart) + call atmos_model_restart(atm_int_state%Atm, timestamp) + + call wrt_atmres_timestamp(atm_int_state,timestamp) + endif endif endif ! @@ -847,20 +885,21 @@ subroutine fcst_finalize(fcst_comp, importState, exportState,clock,rc) 'final time does not match expected ending time', WARNING) !*** write restart file - - call get_date (atm_int_state%Time_atmos, date(1), date(2), date(3), & + if( restart_endfcst ) then + call get_date (atm_int_state%Time_atmos, date(1), date(2), date(3), & date(4), date(5), date(6)) - call mpp_open( unit, 'RESTART/coupler.res', nohdrs=.TRUE. ) - if (mpp_pe() == mpp_root_pe())then - write( unit, '(i6,8x,a)' )calendar_type, & + call mpp_open( unit, 'RESTART/coupler.res', nohdrs=.TRUE. ) + if (mpp_pe() == mpp_root_pe())then + write( unit, '(i6,8x,a)' )calendar_type, & '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - write( unit, '(6i6,8x,a)' )date_init, & + write( unit, '(6i6,8x,a)' )date_init, & 'Model start time: year, month, day, hour, minute, second' - write( unit, '(6i6,8x,a)' )date, & + write( unit, '(6i6,8x,a)' )date, & 'Current model time: year, month, day, hour, minute, second' + endif + call mpp_close(unit) endif - call mpp_close(unit) ! call diag_manager_end(atm_int_state%Time_atmos ) diff --git a/module_fv3_config.F90 b/module_fv3_config.F90 index eb1bb2036..2e0fd8883 100644 --- a/module_fv3_config.F90 +++ b/module_fv3_config.F90 @@ -13,11 +13,10 @@ module module_fv3_config implicit none ! - - integer :: restart_interval -! integer :: nfhout, nfhout_hf, nsout, dt_atmos integer :: nfhmax, nfhmax_hf, first_kdt + integer :: num_restart_interval + integer :: frestart(1000) type(ESMF_Alarm) :: alarm_output_hf, alarm_output type(ESMF_TimeInterval) :: output_hfmax type(ESMF_TimeInterval) :: output_interval,output_interval_hf @@ -25,7 +24,9 @@ module module_fv3_config logical :: cpl, cplprint_flag logical :: quilting, output_1st_tstep_rst logical :: force_date_from_configure + logical :: restart_endfcst ! + integer,dimension(:),allocatable :: restart_interval character(esmf_maxstr),dimension(:),allocatable :: filename_base character(17) :: calendar=' ' integer :: calendar_type = -99 From 63ddb607a967758f7ea75b0d860541e58c1f6925 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 6 Jan 2020 16:07:42 -0500 Subject: [PATCH 06/21] regression test for iccn=1 and iccn=2 --- ccpp/physics | 2 +- gfsphysics/physics/m_micro_driver.F90 | 32 ++++++++++----------------- gfsphysics/physics/micro_mg2_0.F90 | 11 ++++----- gfsphysics/physics/micro_mg3_0.F90 | 11 ++++----- 4 files changed, 25 insertions(+), 31 deletions(-) diff --git a/ccpp/physics b/ccpp/physics index b30d0ce83..e9e685055 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit b30d0ce83f25e80616770dd433451b99eccb3a57 +Subproject commit e9e685055dd95ddcac721e1eaf4a878658e12749 diff --git a/gfsphysics/physics/m_micro_driver.F90 b/gfsphysics/physics/m_micro_driver.F90 index 9d4d3a318..203fe240e 100644 --- a/gfsphysics/physics/m_micro_driver.F90 +++ b/gfsphysics/physics/m_micro_driver.F90 @@ -14,7 +14,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & &, CLLS_io, KCBL & &, CLDREFFL, CLDREFFI, CLDREFFR, CLDREFFS & &, CLDREFFG, aerfld_i & - &, aero_in, naai_i, npccn_i, iccn & + &, naai_i, npccn_i, iccn & &, skip_macro & ! &, skip_macro, cn_prc2, cn_snr & &, lprnt, alf_fac, qc_min, pdfflag & @@ -62,7 +62,8 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & integer, parameter :: ncolmicro = 1 integer,intent(in) :: im, ix,lm, ipr, kdt, fprcp, pdfflag - logical,intent(in) :: flipv, aero_in, skip_macro, lprnt, iccn + logical,intent(in) :: flipv, skip_macro, lprnt + integer,intent(in) :: iccn real (kind=kind_phys), intent(in):: dt_i, alf_fac, qc_min(2) real (kind=kind_phys), dimension(ix,lm),intent(in) :: & @@ -304,7 +305,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,ll) radheat(i,k) = lwheat_i(i,ll) + swheat_i(i,ll) rhc(i,k) = rhc_i(i,ll) - if (iccn) then + if (iccn == 1) then CDNC_NUC(i,k) = npccn_i(i,ll) INC_NUC(i,k) = naai_i (i,ll) endif @@ -365,7 +366,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & temp(i,k) = t_io(i,k) radheat(i,k) = lwheat_i(i,k) + swheat_i(i,k) rhc(i,k) = rhc_i(i,k) - if (iccn) then + if (iccn == 1) then CDNC_NUC(i,k) = npccn_i(i,k) INC_NUC(i,k) = naai_i (i,k) endif @@ -528,7 +529,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & NCPL(i,l) = MAX( NCPL(i,l), 0.) NCPI(i,l) = MAX( NCPI(i,l), 0.) RAD_CF(i,l) = max(0.0, min(CLLS(i,l)+CLCN(i,l), 1.0)) - if (.not. iccn) then + if (iccn.ne.1) then CDNC_NUC(i,l) = 0.0 INC_NUC(i,l) = 0.0 endif @@ -583,7 +584,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! allocate(AERMASSMIX(IM,LM,15)) - if ( aero_in ) then + if (iccn == 2) then AERMASSMIX(:,:,1:ntrcaer) = aerfld_i(:,:,1:ntrcaer) else AERMASSMIX(:,:,1:5) = 1.e-6 @@ -766,12 +767,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & tauxr8 = ter8(K) endif -! if(aero_in) then AeroAux = AeroProps(I, K) -! else -! call init_Aer(AeroAux) -! call init_Aer(AeroAux_b) -! endif pfrz_inc_r8(k) = 0.0 rh1_r8 = 0.0 !related to cnv_dql_dt, needed to changed soon @@ -846,7 +842,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & SC_ICE(i,k) = ((tice-temp(i,k))*tx1 + (temp(i,k)-t_ice_all)*rhc(i,k)) & * t_ice_denom endif - if (.not. iccn) then + if (iccn.ne.1) then CDNC_NUC(I,k) = npccninr8(k) INC_NUC (I,k) = naair8(k) endif @@ -981,7 +977,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! temp(i,k) = th1(i,k) * PK(i,k) RAD_CF(i,k) = min(CLLS(i,k)+CLCN(i,k), 1.0) ! - if (.not. iccn) then + if (iccn.ne.1) then if (PFRZ(i,k) > 0.0) then INC_NUC(i,k) = INC_NUC(i,k) * PFRZ(i,k) NHET_NUC(i,k) = NHET_NUC(i,k) * PFRZ(i,k) @@ -1130,11 +1126,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & endif -! if(aero_in) then - AeroAux = AeroProps(I, K) -! else -! call init_Aer(AeroAux) -! end if + AeroAux = AeroProps(I, K) call getINsubset(1, AeroAux, AeroAux_b) naux = AeroAux_b%nmods if (nbincontactdust < naux) then @@ -1348,7 +1340,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & drout2, dsout2, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, & & lev_sed_strt) ! LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) @@ -1482,7 +1474,7 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & & qgout2, ngout2, dgout2, freqg, & & freqs, freqr, & & nfice, qcrat, & - & prer_evap, xlat(i), xlon(i), lprint, iccn, aero_in, & + & prer_evap, xlat(i), xlon(i), lprint, iccn, & & lev_sed_strt) LS_PRC2(I) = max(1000.*(prectr8(1)-precir8(1)), 0.0) diff --git a/gfsphysics/physics/micro_mg2_0.F90 b/gfsphysics/physics/micro_mg2_0.F90 index 325a2dbbe..82d34ed9f 100644 --- a/gfsphysics/physics/micro_mg2_0.F90 +++ b/gfsphysics/physics/micro_mg2_0.F90 @@ -403,7 +403,7 @@ subroutine micro_mg_tend ( & drout2, dsout2, & freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) + prer_evap, xlat, xlon, lprnt, iccn, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -474,7 +474,8 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt, iccn, aero_in + logical, intent(in) :: lprnt + integer, intent(in) :: iccn ! used for scavenging @@ -1123,7 +1124,7 @@ subroutine micro_mg_tend ( & enddo enddo - if(iccn) then + if(iccn == 1) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1162,7 +1163,7 @@ subroutine micro_mg_tend ( & ncal = zero end where - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1177,7 +1178,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (aero_in) then + elseif (iccn == 2) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then diff --git a/gfsphysics/physics/micro_mg3_0.F90 b/gfsphysics/physics/micro_mg3_0.F90 index cbd25370a..db0870bee 100644 --- a/gfsphysics/physics/micro_mg3_0.F90 +++ b/gfsphysics/physics/micro_mg3_0.F90 @@ -505,7 +505,7 @@ subroutine micro_mg_tend ( & !--ag freqs, freqr, & nfice, qcrat, & - prer_evap, xlat, xlon, lprnt, iccn, aero_in, nlball) + prer_evap, xlat, xlon, lprnt, iccn, nlball) ! Constituent properties. use micro_mg_utils, only: mg_liq_props, & @@ -593,7 +593,8 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: liqcldf(mgncol,nlev) ! liquid cloud fraction (no units) real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) - logical, intent(in) :: lprnt, iccn, aero_in + logical, intent(in) :: lprnt + integer, intent(in) :: iccn ! used for scavenging @@ -1441,7 +1442,7 @@ subroutine micro_mg_tend ( & enddo enddo ! - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol npccn(i,k) = npccnin(i,k) @@ -1495,7 +1496,7 @@ subroutine micro_mg_tend ( & enddo enddo - if (iccn) then + if (iccn == 1) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then @@ -1510,7 +1511,7 @@ subroutine micro_mg_tend ( & endif enddo enddo - elseif (aero_in) then + elseif (iccn == 2) then do k=1,nlev do i=1,mgncol if (t(i,k) < icenuct) then From 58dd381f4d69fc48756b6891d0cb3979da9533f9 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 6 Jan 2020 16:10:39 -0500 Subject: [PATCH 07/21] regression test for iccn=1 and iccn=2 --- ccpp/physics | 2 +- ccpp/suites/suite_FV3_GFS_v16_csawmg.xml | 94 ++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 1 deletion(-) create mode 100644 ccpp/suites/suite_FV3_GFS_v16_csawmg.xml diff --git a/ccpp/physics b/ccpp/physics index e9e685055..47ecb07ac 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit e9e685055dd95ddcac721e1eaf4a878658e12749 +Subproject commit 47ecb07ac0e7dfee9537fa5a013b7bf1e9ed9b7c diff --git a/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml new file mode 100644 index 000000000..147d00891 --- /dev/null +++ b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml @@ -0,0 +1,94 @@ + + + + + + + GFS_time_vary_pre + GFS_rrtmg_setup + GFS_rad_time_vary + GFS_phys_time_vary + + + + + GFS_suite_interstitial_rad_reset + GFS_rrtmg_pre + rrtmg_sw_pre + rrtmg_sw + rrtmg_sw_post + rrtmg_lw_pre + rrtmg_lw + rrtmg_lw_post + GFS_rrtmg_post + + + + + GFS_suite_interstitial_phys_reset + GFS_suite_stateout_reset + get_prs_fv3 + GFS_suite_interstitial_1 + GFS_surface_generic_pre + GFS_surface_composites_pre + dcyc2t3 + GFS_surface_composites_inter + GFS_suite_interstitial_2 + + + + sfc_diff + GFS_surface_loop_control_part1 + sfc_nst_pre + sfc_nst + sfc_nst_post + lsm_noah + sfc_sice + GFS_surface_loop_control_part2 + + + + GFS_surface_composites_post + dcyc2t3_post + sfc_diag + sfc_diag_post + GFS_surface_generic_post + GFS_PBL_generic_pre + satmedmfvdif + GFS_PBL_generic_post + GFS_GWD_generic_pre + cires_ugwp + cires_ugwp_post + rayleigh_damp + GFS_suite_stateout_update + ozphys_2015 + h2ophys + GFS_DCNV_generic_pre + get_phi_fv3 + GFS_suite_interstitial_3 + cs_conv_pre + cs_conv + cs_conv_post + GFS_DCNV_generic_post + GFS_SCNV_generic_pre + samfshalcnv + samfshalcnv_post + GFS_SCNV_generic_post + GFS_suite_interstitial_4 + cnvc90 + GFS_MP_generic_pre + m_micro_pre + m_micro + m_micro_post + cs_conv_aw_adj + GFS_MP_generic_post + maximum_hourly_diagnostics + + + + + GFS_stochastics + + + + From ba967ffe67c563301c5ac00e3a5d385527fd6ede Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Fri, 17 Jan 2020 14:16:03 -0500 Subject: [PATCH 08/21] MERRA2 consistent radiation pass regression tests --- gfsphysics/GFS_layer/GFS_diagnostics.F90 | 79 +- gfsphysics/GFS_layer/GFS_radiation_driver.F90 | 21 +- gfsphysics/physics/aerclm_def.f | 16 +- gfsphysics/physics/aerinterp.f90 | 320 +- gfsphysics/physics/radiation_aerosols.f | 2867 +++++------------ 5 files changed, 1108 insertions(+), 2195 deletions(-) diff --git a/gfsphysics/GFS_layer/GFS_diagnostics.F90 b/gfsphysics/GFS_layer/GFS_diagnostics.F90 index ec6c1233b..80d2adc3f 100644 --- a/gfsphysics/GFS_layer/GFS_diagnostics.F90 +++ b/gfsphysics/GFS_layer/GFS_diagnostics.F90 @@ -700,8 +700,85 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,16) ExtDiag(idx)%data(nb)%var21 => IntDiag(nb)%fluxr(:,7) enddo -! if(mpp_pe()==mpp_root_pe())print *,'in gfdl_diag_register,af TEMP_avelct,idx=',idx +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'AOD_550' + ExtDiag(idx)%desc = 'total aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,34) + enddo + +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'DU_AOD_550' + ExtDiag(idx)%desc = 'dust aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,35) + enddo + +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'BC_AOD_550' + ExtDiag(idx)%desc = 'soot aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,36) + enddo + +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'OC_AOD_550' + ExtDiag(idx)%desc = 'waso aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,37) + enddo +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'SU_AOD_550' + ExtDiag(idx)%desc = 'suso aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,38) + enddo + +!--- aerosol diagnostics --- + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'SS_AOD_550' + ExtDiag(idx)%desc = 'salt aerosol optical depth at 550 nm' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%intpl_method = 'bilinear' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%fluxr(:,39) + enddo + +! ! !--- accumulated diagnostics --- do num = 1,NFXR diff --git a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 index 5215f1b2b..7200d4616 100644 --- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 @@ -1561,7 +1561,8 @@ subroutine GFS_radiation_driver & !check print *,' in grrad : calling setaer ' call setaer (plvl, plyr, prslk1, tvly, rhly, Sfcprop%slmsk, & ! --- inputs - tracer1, Grid%xlon, Grid%xlat, IM, LMK, LMP, & + tracer1, Tbd%aer_nm, & + Grid%xlon, Grid%xlat, IM, LMK, LMP, & Model%lsswr,Model%lslwr, & faersw,faerlw,aerodp) ! --- outputs @@ -2039,12 +2040,18 @@ subroutine GFS_radiation_driver & if (Model%lssav) then if (Model%lsswr) then do i=1,im - Diag%fluxr(i,34) = Diag%fluxr(i,34) + Model%fhswr*aerodp(i,1) ! total aod at 550nm - Diag%fluxr(i,35) = Diag%fluxr(i,35) + Model%fhswr*aerodp(i,2) ! DU aod at 550nm - Diag%fluxr(i,36) = Diag%fluxr(i,36) + Model%fhswr*aerodp(i,3) ! BC aod at 550nm - Diag%fluxr(i,37) = Diag%fluxr(i,37) + Model%fhswr*aerodp(i,4) ! OC aod at 550nm - Diag%fluxr(i,38) = Diag%fluxr(i,38) + Model%fhswr*aerodp(i,5) ! SU aod at 550nm - Diag%fluxr(i,39) = Diag%fluxr(i,39) + Model%fhswr*aerodp(i,6) ! SS aod at 550nm +! Diag%fluxr(i,34) = Diag%fluxr(i,34) + Model%fhswr*aerodp(i,1) ! total aod at 550nm +! Diag%fluxr(i,35) = Diag%fluxr(i,35) + Model%fhswr*aerodp(i,2) ! DU aod at 550nm +! Diag%fluxr(i,36) = Diag%fluxr(i,36) + Model%fhswr*aerodp(i,3) ! BC aod at 550nm +! Diag%fluxr(i,37) = Diag%fluxr(i,37) + Model%fhswr*aerodp(i,4) ! OC aod at 550nm +! Diag%fluxr(i,38) = Diag%fluxr(i,38) + Model%fhswr*aerodp(i,5) ! SU aod at 550nm +! Diag%fluxr(i,39) = Diag%fluxr(i,39) + Model%fhswr*aerodp(i,6) ! SS aod at 550nm + Diag%fluxr(i,34) = aerodp(i,1) ! total aod at 550nm + Diag%fluxr(i,35) = aerodp(i,2) ! DU aod at 550nm + Diag%fluxr(i,36) = aerodp(i,3) ! BC aod at 550nm + Diag%fluxr(i,37) = aerodp(i,4) ! OC aod at 550nm + Diag%fluxr(i,38) = aerodp(i,5) ! SU aod at 550nm + Diag%fluxr(i,39) = aerodp(i,6) ! SS aod at 550nm enddo endif diff --git a/gfsphysics/physics/aerclm_def.f b/gfsphysics/physics/aerclm_def.f index 6729237d8..84852a1de 100644 --- a/gfsphysics/physics/aerclm_def.f +++ b/gfsphysics/physics/aerclm_def.f @@ -2,22 +2,22 @@ module aerclm_def use machine , only : kind_phys implicit none -! only read monthly merra2 data for m-1, m, m+1 - integer, parameter :: levsaer=45, latsaer=91, lonsaer=144 - integer, parameter :: lmerra=72, ntrcaerm=15, timeaer=12 + integer, parameter :: levsaer=50, ntrcaerm=15, timeaer=12 + integer :: latsaer, lonsaer, ntrcaer - integer :: ntrcaer character*10 :: specname(ntrcaerm) - real (kind=kind_phys):: aer_lat(latsaer), aer_lon(lonsaer) - & ,aer_time(13) - real (kind=4), allocatable, dimension(:,:,:,:,:) :: aerin + 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_pres + real (kind=kind_phys), allocatable, dimension(:,:,:,:,:) :: aerin data aer_time/15.5, 45., 74.5, 105., 135.5, 166., 196.5, & 227.5, 258., 288.5, 319., 349.5, 380.5/ data specname /'DU001','DU002','DU003','DU004','DU005', & 'SS001','SS002','SS003','SS004','SS005','SO4', - & 'BCPHOBIC','BCPHILIC','OCPHILIC','OCPHOBIC'/ + & 'BCPHOBIC','BCPHILIC','OCPHOBIC','OCPHILIC'/ end module aerclm_def diff --git a/gfsphysics/physics/aerinterp.f90 b/gfsphysics/physics/aerinterp.f90 index 8d5603b83..8d35efbcf 100644 --- a/gfsphysics/physics/aerinterp.f90 +++ b/gfsphysics/physics/aerinterp.f90 @@ -1,165 +1,181 @@ SUBROUTINE read_aerdata (me, master, iflip, idate ) - - use machine, only: kind_phys + use machine, only: kind_phys, kind_io4, kind_io8 use aerclm_def use netcdf !--- in/out integer, intent(in) :: me, master, iflip, idate(4) - !--- locals - integer :: ncid, varid - integer :: i, j, k, n, ii, ijk, imon, klev - character :: fname*50, mn*2, fldname*10 + integer :: ncid, varid, ndims, dim1, dim2, dim3, hmx + integer :: i, j, k, n, ii, imon, klev + character :: fname*50, mn*2, vname*10 logical :: file_exist - real(kind=4), allocatable, dimension(:,:,:) :: ps_clm - real(kind=4), allocatable, dimension(:,:,:,:) :: delp_clm - real(kind=4), allocatable, dimension(:,:,:,:) :: aer_clm - real(kind=4), allocatable, dimension(:,:,:,:) :: airden_clm - real(kind=4), allocatable, dimension(:) :: pres_tmp - - allocate (delp_clm(lonsaer,latsaer,lmerra,1)) - allocate (aer_clm(lonsaer,latsaer,lmerra,1)) - allocate (airden_clm(lonsaer,latsaer,lmerra,1)) - allocate (ps_clm(lonsaer,latsaer,1)) - allocate (pres_tmp(lmerra)) - -! allocate aerclm_def arrays: aerin and aer_pres - allocate (aerin(lonsaer,latsaer,levsaer,ntrcaer,timeaer)) - allocate (aer_pres(lonsaer,latsaer,levsaer,timeaer)) + integer, allocatable :: invardims(:) + real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff + real(kind=kind_io4),allocatable,dimension(:,:,:,:):: buffx + real(kind=kind_io4),allocatable,dimension(:,:) :: pres_tmp + real(kind=kind_io8),allocatable,dimension(:) :: aer_lati + real(kind=kind_io8),allocatable,dimension(:) :: aer_loni +! +!! =================================================================== if (me == master) then if ( iflip == 0 ) then ! data from toa to sfc - print *, "EJ, GFS is top-down" + print *, "GFS is top-down" else - print *, "EJ, GFS is bottom-up" + print *, "GFS is bottom-up" endif endif +! +!! =================================================================== +!! fetch dim spec and lat/lon from m01 data set +!! =================================================================== + fname=trim("aeroclim.m"//'01'//".nc") + inquire (file = fname, exist = file_exist) + if (.not. file_exist ) then + print *, 'fname not found, abort' + stop 8888 + endif + call nf_open(fname , NF90_NOWRITE, ncid) + + vname = trim(specname(1)) + call nf_inq_varid(ncid, vname, varid) + call nf_inq_varndims(ncid, varid, ndims) + + if(.not. allocated(invardims)) allocate(invardims(3)) + call nf_inq_vardimid(ncid,varid,invardims) + call nf_inq_dimlen(ncid, invardims(1), dim1) + call nf_inq_dimlen(ncid, invardims(2), dim2) + call nf_inq_dimlen(ncid, invardims(3), dim3) + +! specify latsaer, lonsaer, hmx + lonsaer = dim1 + latsaer = dim2 + hmx = int(dim1/2) ! to swap long from W-E to E-W + + if(me==master) then + print *, 'MERRA2 dim: ',dim1, dim2, dim3 + endif + +! allocate arrays + if (.not. allocated(aer_loni)) then + allocate (aer_loni(lonsaer)) + allocate (aer_lati(latsaer)) + endif + + if (.not. allocated(aer_lat)) then + allocate(aer_lat(latsaer)) + allocate(aer_lon(lonsaer)) + allocate(aerin(lonsaer,latsaer,levsaer,ntrcaerm,timeaer)) + allocate(aer_pres(lonsaer,latsaer,levsaer,timeaer)) + endif +! construct lat/lon array + call nf_inq_varid(ncid, 'lat', varid) + call nf_get_var(ncid, varid, aer_lati) + call nf_inq_varid(ncid, 'lon', varid) + call nf_get_var(ncid, varid, aer_loni) + + do i = 1, hmx ! flip from (-180,180) to (0,360) + if(aer_loni(i)<0.) aer_loni(i)=aer_loni(i)+360. + aer_lon(i+hmx) = aer_loni(i) + aer_lon(i) = aer_loni(i+hmx) + enddo + + do i = 1, latsaer + aer_lat(i) = aer_lati(i) + enddo + + call nf_close(ncid) + +! allocate local working arrays + if (.not. allocated(buff)) then + allocate (buff(lonsaer, latsaer, dim3)) + allocate (pres_tmp(lonsaer,dim3)) + endif + if (.not. allocated(buffx)) then + allocate (buffx(lonsaer, latsaer, dim3,1)) + endif + +!! =================================================================== +!! loop thru m01 - m12 for aer/pres array +!! =================================================================== do imon = 1, timeaer - !ijk = imon + idate(2)+int(idate(3)/16)-2 - !if ( ijk .le. 0 ) ijk = 12 - !if ( ijk .eq. 13 ) ijk = 1 - !if ( ijk .eq. 14 ) ijk = 2 write(mn,'(i2.2)') imon - fname=trim("merra2C.aerclim.2003-2014.m"//mn//".nc") - if (me == master) print *, "EJ,aerosol climo:", fname, & + fname=trim("aeroclim.m"//mn//".nc") + if (me == master) print *, "aerosol climo:", fname, & "for imon:",imon,idate inquire (file = fname, exist = file_exist) if ( file_exist ) then if (me == master) print *, & - "EJ, aerosol climo found; proceed the run" + "aerosol climo found; proceed the run" else - print *,"EJ, Error! aerosol climo not found; abort the run" + print *,"Error! aerosol climo not found; abort the run" stop 555 endif - call nf_open(fname, nf_NOWRITE, ncid) - -! merra2 data is top down -! for GFS, iflip 0: toa to sfc; 1: sfc to toa - -! read aerosol mixing ratio arrays (kg/kg) -! construct 4-d aerosol mass concentration (kg/m3) - call nf_inq_varid(ncid, 'AIRDENS', varid) - call nf_get_var(ncid, varid, airden_clm) -! if(me==master) print *, "EJ, read airdens", airden_clm(1,1,:,1) + call nf_open(fname , nf90_NOWRITE, ncid) - do ii = 1, ntrcaer - fldname=specname(ii) - call nf_inq_varid(ncid, fldname, varid) - call nf_get_var(ncid, varid, aer_clm) -! if(me==master) print *, "EJ, read ", fldname, aer_clm(1,1,:,1) - do i = 1, lonsaer - do j = 1, latsaer - 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 = ( lmerra - k ) + 1 - endif - aerin(i,j,k,ii,imon) = aer_clm(i,j,klev,1)*airden_clm(i,j,klev,1) - enddo !k-loop (lev) - enddo !j-loop (lat) - enddo !i-loop (lon) - enddo !ii-loop (ntrac) - -! aer_clm is top-down (following MERRA2) -! aerin is bottom-up (following GFS) - -! if ( imon == 1 .and. me == master ) then -! print *, 'EJ, du1(1,1) :', aerin(1,1,:,1,imon) -! endif - -! construct 3-d pressure array (Pa) - call nf_inq_varid(ncid, "PS", varid) - call nf_get_var(ncid, varid, ps_clm) +! ====> construct 3-d pressure array (Pa) call nf_inq_varid(ncid, "DELP", varid) - call nf_get_var(ncid, varid, delp_clm) - -! if ( imon == 1 .and. me == master ) then -! print *, 'EJ, ps_clm:', ps_clm(1,1,1) -! print *, 'EJ, delp_clm:', delp_clm(1,1,:,1) -! endif + call nf_get_var(ncid, varid, buff) - do i = 1, lonsaer do j = 1, latsaer + do i = 1, lonsaer +! constract pres_tmp (top-down), note input is top-down + pres_tmp(i,1) = 0. + do k=2, dim3 + pres_tmp(i,k) = pres_tmp(i,k-1)+buff(i,j,k) + enddo !k-loop + enddo !i-loop (lon) -! constract pres_tmp (top-down) - pres_tmp(1) = 0. - do k=2, lmerra - pres_tmp(k) = pres_tmp(k-1) + delp_clm(i,j,k,1) - enddo -! if (imon==1 .and. me==master .and. i==1 .and. j==1 ) then -! print *, 'EJ, pres_tmp:', pres_tmp(:) -! endif - -! extract pres_tmp to fill aer_pres +! 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 = ( lmerra - k ) + 1 + klev = ( dim3 - k ) + 1 endif - aer_pres(i,j,k,imon)= pres_tmp(klev) + do i = 1, hmx + aer_pres(i+hmx,j,k,imon)= 1.d0*pres_tmp(i,klev) + aer_pres(i,j,k,imon) = 1.d0*pres_tmp(i+hmx,klev) + enddo !i-loop (lon) enddo !k-loop (lev) -! if (imon==1 .and. me==master .and. i==1 .and. j==1 ) then -! print *, 'EJ, aer_pres:', aer_pres(i,j,:,imon) -! endif - enddo !j-loop (lat) - enddo !i-loop (lon) -! if (imon==1 .and. me==master ) then -! print *, 'EJx, aer_pres_i1:',(aer_pres(1,1:180,levsaer,imon) ) -! endif +! ====> 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)) + call nf_inq_varid(ncid, vname, varid) + call nf_get_var(ncid, varid, buffx) -! construct lat/lon array - if (imon == 1 ) then - call nf_inq_varid(ncid, "lat", varid) - call nf_get_var(ncid, varid, aer_lat) - call nf_inq_varid(ncid, "lon", varid) - call nf_get_var(ncid, varid, aer_lon) - do i = 1, lonsaer - if(aer_lon(i) < 0.) aer_lon(i) = aer_lon(i) + 360. - enddo -! if (imon==1 .and. me == master) then -! print *, "EJ, lat:", aer_lat(:) -! print *, "EJ, lon:", aer_lon(:) -! endif - endif + do j = 1, latsaer + 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 = ( dim3 - k ) + 1 + endif + do i = 1, hmx + aerin(i+hmx,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1) + aerin(i,j,k,ii,imon) = 1.d0*buffx(i+hmx,j,klev,1) + enddo !i-loop (lon) + enddo !k-loop (lev) + enddo !j-loop (lat) + + ENDDO ! ii-loop (ntracaerm) ! close the file call nf_close(ncid) enddo !imon-loop - !--- - deallocate (ps_clm, delp_clm, pres_tmp, aer_clm, airden_clm ) - if (me == master) then - write(*,*) 'Reading in GOCART aerosols data' - endif + deallocate (aer_loni, aer_lati) + deallocate (buff, pres_tmp) + deallocate (buffx) END SUBROUTINE read_aerdata ! @@ -197,11 +213,6 @@ SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, & ddy(j) = 1.0 endif -! if (me == master .and. j<= 3) then -! print *,'EJj,',j,' dlat=',dlat(j),' jindx12=',jindx1(j),& -! jindx2(j),' aer_lat=',aer_lat(jindx1(j)), & -! aer_lat(jindx2(j)),' ddy=',ddy(j) -! endif ENDDO DO J=1,npts @@ -220,11 +231,6 @@ SUBROUTINE setindxaer(npts,dlat,jindx1,jindx2,ddy,dlon, & else ddx(j) = 1.0 endif -! if (me == master .and. j<= 3) then -! print *,'EJi,',j,' dlon=',dlon(j),' iindx12=',iindx1(j),& -! iindx2(j),' aer_lon=',aer_lon(iindx1(j)), & -! aer_lon(iindx2(j)),' ddx=',ddx(j) -! endif ENDDO RETURN @@ -248,7 +254,8 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & integer IDAT(8),JDAT(8) ! real(kind=kind_phys) DDY(npts), ddx(npts),ttt - real(kind=kind_phys) aerout(npts,lev,ntrcaer),aerpm(npts,levsaer,ntrcaer) + 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) RINC(5), rjday integer jdow, jdoy, jday @@ -269,7 +276,6 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & else CALL W3MOVDAT(RINC,IDAT,JDAT) endif -! if(me==master) print *,'EJ, IDAT ',IDAT(1:3), IDAT(5) ! jdow = 0 jdoy = 0 @@ -290,15 +296,8 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & tx1 = (aer_time(n2) - rjday) / (aer_time(n2) - aer_time(n1)) tx2 = 1.0 - tx1 if (n2 > 12) n2 = n2 -12 -! if(me==master)print *,'EJ,rjday=',rjday, ';aer_time,tx1,tx=' & -! , aer_time(n1),aer_time(n2),tx1,tx2,n1,n2 -! -! if(me==master) then -! DO L=1,levsaer -! print *,'EJ,aerin(n1,n2)=',L,aerin(1,1,L,1,n1),aerin(1,1,L,1,n2) -! ENDDO -! endif +! DO L=1,levsaer DO J=1,npts J1 = JINDX1(J) @@ -321,49 +320,38 @@ SUBROUTINE aerinterpol(me,master,npts,IDATE,FHOUR,jindx1,jindx2, & +tx2*(TEMI*TEMJ*aer_pres(I1,J1,L,n2)+DDX(j)*DDY(J)*aer_pres(I2,J2,L,n2) & +TEMI*DDY(j)*aer_pres(I1,J2,L,n2)+DDX(j)*TEMJ*aer_pres(I2,J1,L,n2)) -! IF(me==master .and. j==1) THEN -! print *, 'EJ,aer/ps:',L,aerpm(j,L,1),aerpres(j,L) -! if(L==1) then -! print *, 'EJ, wgt:',TEMI*TEMJ,DDX(j)*DDY(J),TEMI*DDY(j),DDX(j)*TEMJ -! print *, 'EJ, aerx:',aerin(I1,J1,L,ii,n1), & -! aerin(I2,J2,L,ii,n1), aerin(I1,J2,L,ii,n1), aerin(I2,J1,L,ii,n1) -! print *, 'EJ, aery:',aerin(I1,J1,L,ii,n2), & -! aerin(I2,J2,L,ii,n2), aerin(I1,J2,L,ii,n2), aerin(I2,J1,L,ii,n2) -! endif -! ENDIF ENDDO ENDDO -! note: input is set to be same as GFS +! don't flip, input is the same direction as GFS (bottom-up) DO J=1,npts DO L=1,lev - if(prsl(j,l).ge.aerpres(j,levsaer)) then + if(prsl(j,L).ge.aerpres(j,1)) then DO ii=1, ntrcaer - aerout(j,l,ii)=aerpm(j,levsaer,ii) + aerout(j,L,ii)=aerpm(j,1,ii) !! sfc level ENDDO - else if(prsl(j,l).le.aerpres(j,1)) then + else if(prsl(j,L).le.aerpres(j,levsaer)) then DO ii=1, ntrcaer - aerout(j,l,ii)=aerpm(j,1,ii) + aerout(j,L,ii)=aerpm(j,levsaer,ii) !! toa top ENDDO else - DO k=levsaer-1,1,-1 - IF(prsl(j,l)>aerpres(j,k)) then + DO k=1, levsaer-1 !! from sfc to toa + IF(prsl(j,L)aerpres(j,k+1)) then i1=k i2=min(k+1,levsaer) exit - end if - end do + ENDIF + ENDDO + temi = prsl(j,L)-aerpres(j,i2) + temj = aerpres(j,i1) - prsl(j,L) + tx1 = temi/(aerpres(j,i1) - aerpres(j,i2)) + tx2 = temj/(aerpres(j,i1) - aerpres(j,i2)) DO ii = 1, ntrcaer - aerout(j,l,ii)=aerpm(j,i1,ii)+(aerpm(j,i2,ii)-aerpm(j,i1,ii))& - /(aerpres(j,i2)-aerpres(j,i1))*(prsl(j,l)-aerpres(j,i1)) -! IF(me==master .and. j==1 .and. ii==1) then -! print *, 'EJ, aerout:',aerout(j,l,ii), aerpm(j,i1,ii), & -! aerpm(j,i2,ii), aerpres(j,i2), aerpres(j,i1), prsl(j,l) -! ENDIF + aerout(j,L,ii)= aerpm(j,i1,ii)*tx1 + aerpm(j,i2,ii)*tx2 ENDDO - endif - ENDDO - ENDDO + endif + ENDDO !L-loop + ENDDO !J-loop ! - RETURN + RETURN END diff --git a/gfsphysics/physics/radiation_aerosols.f b/gfsphysics/physics/radiation_aerosols.f index 37364b8be..339b991f0 100644 --- a/gfsphysics/physics/radiation_aerosols.f +++ b/gfsphysics/physics/radiation_aerosols.f @@ -25,11 +25,10 @@ ! ! ! 'setaer' -- mapping aeros profile, compute aeros opticals ! ! inputs: ! -! (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat, ! +! (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat, ! ! IMAX,NLAY,NLP1, lsswr,lslwr, ! ! outputs: ! -! (aerosw,aerolw,tau_gocart) ! -!! (aerosw,aerolw,aerodp) ! +! (aerosw,aerolw,aerodp) ! ! ! ! ! ! external modules referenced: ! @@ -100,6 +99,9 @@ ! jun 2018 --- h-m lin and y-t hou updated spectral band ! ! mapping method for aerosol optical properties. controled by ! ! internal variable lmap_new through namelist variable iaer. ! +! may 2019 --- sarah lu, restore the gocart option, allowing ! +! aerosol ext, ssa, asy determined from MERRA2 monthly climo ! +! with new spectral band mapping method ! ! ! ! references for opac climatological aerosols: ! ! hou et al. 2002 (ncep office note 441) ! @@ -107,6 +109,11 @@ ! ! ! references for gocart interactive aerosols: ! ! chin et al., 2000 - jgr, v105, 24671-24687 ! +! colarco et al., 2010 - jgr, v115, D14207 ! +! ! +! references for merra2 aerosol reanalysis: ! +! randles et al., 2017 - jclim, v30, 6823-6850 ! +! buchard et al., 2017 - jclim, v30, 6851-6871 ! ! ! ! references for stratosperic volcanical aerosols: ! ! sato et al. 1993 - jgr, v98, d12, 22987-22994 ! @@ -139,6 +146,12 @@ !! \cite hess_et_al_1998 !! - GOCART interactive aerosols: !! Chin et al., 2000 \cite chin_et_al_2000 +!! Colarco et al., 2010 - jgr, v115, D14207\cite colarco_et_al_2010 +!! +!! - MERRA2 aerosol reanalysis: +!! Randles et al., 2017 - jclim, v30, 6823-6850\cite randles_et_al_2017 +!! Buchard et al., 2017 - jclim, v30, 6851-6871\cite buchard_et_al_2017 +!! !! - Stratospheric volcanical aerosols: !! Sato et al. 1993 \cite sato_et_al_1993 !========================================! @@ -156,7 +169,8 @@ module module_radiation_aerosols ! use module_radlw_parameters, only : NBDLW, wvnlw1, wvnlw2 ! use funcphys, only : fpkap - use gfs_phy_tracer_config, only : gfs_phy_tracer, trcindx + use aerclm_def, only : ntrcaer + ! implicit none ! @@ -393,24 +407,20 @@ module module_radiation_aerosols ! ! --------------------------------------------------------------------- ! ! section-4 : module variables for gocart aerosol optical properties ! ! --------------------------------------------------------------------- ! - !> \name module variables for gocart aerosol optical properties ! --- parameters and constants: -! - KCM, KCM1, KCM2 are determined from subroutine 'set_aerspc' !> num of bands for aer data (gocart) - integer, parameter :: KAERBND=61 + integer, parameter :: KAERBNDD=61 + integer, parameter :: KAERBNDI=56 !> num of rh levels for rh-dep components integer, parameter :: KRHLEV =36 -!* integer, parameter :: KCM1 = 8 ! num of rh independent aer !species -!* integer, parameter :: KCM2 = 5 ! num of rh dependent aer species -!* integer, parameter :: KCM = KCM1 + KCM2 -!> num of rh indep aerosols (set in subr set_aerspc) - integer, save :: KCM1 = 0 -!> num of rh dep aerosols (set in subr set_aerspc) - integer, save :: KCM2 = 0 -!> =KCM1+KCM2 (set in subr set_aerspc) - integer, save :: KCM +!> num of gocart rh indep aerosols + integer, parameter :: KCM1 = 5 +!> num of gocart rh dep aerosols + integer, parameter :: KCM2 = 10 +!> num of gocart aerosols + integer, parameter :: KCM = KCM1 + KCM2 real (kind=kind_phys), dimension(KRHLEV) :: rhlev_grt & data rhlev_grt (:)/ .00, .05, .10, .15, .20, .25, .30, .35, & @@ -418,252 +428,42 @@ module module_radiation_aerosols ! & .83, .84, .85, .86, .87, .88, .89, .90, .91, .92, .93, & & .94, .95, .96, .97, .98, .99 / -! --- the following arrays are allocate and setup in subr 'gocrt_aerinit' -! ------ gocart aerosol specification ------ -! => transported aerosol species: -! DU (5-bins) -! SS (4 bins for climo mode and 5 bins for fcst mode) -! SU (dms, so2, so4, msa) -! OC (phobic, philic) and BC (phobic, philic) -! => species and lumped species for aerosol optical properties -! DU (5-bins, with 4 sub-groups in the submicron bin ) -! SS (ssam for submicron, sscm for coarse mode) -! SU (so4) -! OC (phobic, philic) and BC (phobic, philic) -! => specification used for aerosol optical properties luts -! DU (8 bins) -! SS (ssam, sscm) -! SU (suso) -! OC (waso) and BC (soot) -! -! - spectral band structure: -! iendwv_grt(KAERBND) - ending wavenumber (cm**-1) for each band -! - relative humidity independent aerosol optical properties: -! ===> species : dust (8 bins) -! rhidext0_grt(KAERBND,KCM1) - extinction coefficient -! rhidssa0_grt(KAERBND,KCM1) - single scattering albedo -! rhidasy0_grt(KAERBND,KCM1) - asymmetry parameter -! - relative humidity dependent aerosol optical properties: -! ===> species : soot, suso, waso, ssam, sscm -! rhdpext0_grt(KAERBND,KRHLEV,KCM2) - extinction coefficient -! rhdpssa0_grt(KAERBND,KRHLEV,KCM2) - single scattering albedo -! rhdpasy0_grt(KAERBND,KRHLEV,KCM2) - asymmetry parameter - -!> spectral band structure: ending wavenumber (\f$cm^-1\f$) for each band - integer, allocatable, dimension(:) :: iendwv_grt -! relative humidity independent aerosol optical properties: -!! species : dust (8 bins) - !> \name relative humidity independent aerosol optical properties: -!! species : dust (8 bins) - -!> extinction coefficient - real (kind=kind_phys),allocatable, dimension(:,:) :: rhidext0_grt -!> single scattering albedo - real (kind=kind_phys),allocatable, dimension(:,:) :: rhidssa0_grt -!> asymmetry parameter - real (kind=kind_phys), allocatable, dimension(:,:) :: rhidasy0_grt +!! species: du001, du002, du003, du004, du005 +! extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw band +! scarhi_grt(KCM1,NSWLWBD) - scattering coefficient for sw+lw band +! ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw band +! asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw band + real (kind=kind_phys),allocatable,save,dimension(:,:) :: & + & extrhi_grt, scarhi_grt, ssarhi_grt, asyrhi_grt ! -! relative humidity dependent aerosol optical properties: -! species : soot, suso, waso, ssam, sscm - !> \name relative humidity dependent aerosol optical properties: -!! species : soot, suso, waso, ssam, sscm +!! species : ss001, ss002, ss003, ss004, ss005, so4, +!! bcphobic, bcphilic, ocphobic, ocphilic +! extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction coefficient for sw+lw band +! scarhd_grt(KRHLEV,KCM2,NSWLWBD) - scattering coefficient for sw+lw band +! ssarhd_grt(KRHLEV,KCM2,NSWLWBD) - single scattering albedo for sw+lw band +! asyrhd_grt(KRHLEV,KCM2,NSWLWBD) - asymmetry parameter for sw+lw band !> extinction coefficient - real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpext0_grt -!> single scattering albedo - real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpssa0_grt -!> asymmetry parameter - real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpasy0_grt - -! - relative humidity independent aerosol optical properties: -! extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band -! ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band -! asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band -! - relative humidity dependent aerosol optical properties: -! extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction coefficient for sw+lw band -! ssarhd_grt(KRHLEV,KCM2,NSWLWBD) - single scattering albedo for sw+lw band -! asyrhd_grt(KRHLEV,KCM2,NSWLWBD) - asymmetry parameter for sw+lw band - -!>\name relative humidity independent aerosol optical properties - -!> extinction coefficient for SW+LW spectral band - real (kind=kind_phys),allocatable,save,dimension(:,:) :: & - & extrhi_grt -!> single scattering albedo for SW+LW spectral band - real (kind=kind_phys),allocatable,save,dimension(:,:) :: & - & ssarhi_grt -!> asymmetry parameter for SW+LW spectral band - real (kind=kind_phys),allocatable,save,dimension(:,:) :: & - & asyrhi_grt - -!> \name relative humidity dependent aerosol optical properties - -!> extinction coefficient for SW+LW spectral band real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: & - & extrhd_grt -!> single scattering albedo for SW+LW band - real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: & - & ssarhd_grt -!> asymmetry parameter for SW+LW band - real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: & - & asyrhd_grt + & extrhd_grt, scarhd_grt, ssarhd_grt, asyrhd_grt -!> \name module variables for gocart aerosol clim data set +!> gocart species + integer, parameter :: num_gc = 5 + character*2 :: gridcomp(num_gc) + integer, dimension (num_gc):: num_radius, radius_lower + integer, dimension (num_gc):: trc_to_aod -! --------------------------------------------------------------------- ! -! section-5 : module variables for gocart aerosol climo data set ! -! --------------------------------------------------------------------- ! -! This version only supports geos3-gocart data set (Jan 2010) -! Modified to support geos4-gocart data set (May 2010) -! -! geos3-gocart vs geos4-gocart -! (1) Use the same module variables -! IMXG,JMXG,KMXG,NMXG,psclmg,dmclmg,geos_rlon,geos_rlat -! (2) Similarity between geos3 and geos 4: -! identical lat/lon grids and aerosol specification; -! direction of vertical index is bottom-up (sfc to toa) -! (3) Difference between geos3 and geos4 -! vertical coordinate (sigma for geos3/hybrid_sigma_pressure for geos4) -! aerosol units (mass concentration for geos3/mixing ratio for geos4) - -!> num of lon-points in geos dataset - integer, parameter :: IMXG = 144 -!> num of lat-points in geos dataset - integer, parameter :: JMXG = 91 -!> num of vertical layers in geos dataset - integer, parameter :: KMXG = 30 -!* integer, parameter :: NMXG = 12 -!> to be determined by set_aerspc - integer, save :: NMXG - - real (kind=kind_phys), parameter :: dltx = 360.0 / float(IMXG) - real (kind=kind_phys), parameter :: dlty = 180.0 / float(JMXG-1) - -! --- the following arrays are allocated and setup in 'rd_gocart_clim' -! - geos-gocart climo data (input dataset) -! psclmg - pressure in cb IMXG*JMXG*KMXG -! dmclmg - aerosol dry mass in g/m3 IMXG*JMXG*KMXG*NMXG -! or aerosol mixing ratio in mol/mol or Kg/Kg - -!> pressure in cb - real (kind=kind_phys),allocatable, save:: psclmg(:,:,:) -!> aerosol dry mass in g/m3 or aerosol mixing ration in mol/mol or Kg/Kg - real (kind=kind_phys),allocatable, save:: dmclmg(:,:,:,:) - -! - geos-gocart lat/lon arrays -!> geos-gocart longitude arrays - real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlon -!> geos-gocart latitude arrays - real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlat - -!> control flag for gocart climo data set: xxxx as default; ver3 for geos3; -!! ver4 for geos4; 0000 for unknown data - character*4, save :: gocart_climo = 'xxxx' - -!> molecular wght of gocart aerosol species - real (kind=kind_io4), allocatable :: molwgt(:) - -! --------------------------------------------------------------------- -! ! -! section-6 : module variables for gocart aerosol scheme options -! ! -! --------------------------------------------------------------------- -! ! - -!> logical parameter for gocart initialization control - logical, save :: lgrtint = .true. - -!> logical parameter for gocart debug print control -! logical, save :: lckprnt = .true. - logical, save :: lckprnt = .false. - -! --- the following index/flag/weight are set up in 'set_aerspc' - -!> merging coefficients for fcst/clim; determined from fdaer - real (kind=kind_phys), save :: ctaer = f_zero ! user specified wgt - -!> option to get fcst gocart aerosol field - logical, save :: get_fcst = .true. -!> option to get clim gocart aerosol field - logical, save :: get_clim = .true. - -! ------ gocart aerosol specification ------ -! => transported aerosol species: -! DU (5-bins) -! SS (4 bins for climo mode and 5 bins for fcst mode) -! SU (dms, so2, so4, msa) -! OC (phobic, philic) and BC (phobic, philic) -! => species and lumped species for aerosol optical properties -! DU (5-bins, with 4 sub-groups in the submicron bin ) -! SS (ssam for submicron, sscm for coarse mode) -! SU (so4) -! OC (phobic, philic) and BC (phobic, philic) -! => specification used for aerosol optical properties luts -! DU (8 bins) -! SS (ssam, sscm) -! SU (suso) -! OC (waso) and BC (soot) -! + data gridcomp /'DU', 'SS', 'SU', 'BC', 'OC'/ + data num_radius /5, 5, 1, 2, 2 / + data radius_lower /1, 6, 11, 12, 14 / + data trc_to_aod /1, 5, 4, 2, 3/ ! dust, soot, waso, suso, ssam -!> index for rh dependent aerosol optical properties (2nd -!! dimension for extrhd_grt, ssarhd_grt, and asyrhd_grt) - integer, save :: isoot, iwaso, isuso, issam, isscm - -! - index for rh independent aerosol optical properties (1st -! dimension for extrhi_grt, ssarhi_grt, and asyrhi_grt) is -! not needed ===> hardwired to 8-bin dust - -! - index for gocart aerosol species to be included in the -! calculations of aerosol optical properties (ext, ssa, asy) -!> index for gocart aerosol species to be included in the -!! calculations of aerosol optical properties (ext, ssa, asy) - type gocart_index_type -! dust - integer :: dust1, dust2, dust3, dust4, dust5 -! sea salt - integer :: ssam, sscm -! sulfate - integer :: suso -! oc - integer :: waso_phobic, waso_philic -! bc - integer :: soot_phobic, soot_philic - endtype - type (gocart_index_type), save :: dm_indx - -!> index for gocart aerosols from prognostic tracer fields - type tracer_index_type -! dust - integer :: du001, du002, du003, du004, du005 -! sea salt - integer :: ss001, ss002, ss003, ss004, ss005 -! sulfate - integer :: so4 -! oc - integer :: ocphobic, ocphilic -! bc - integer :: bcphobic, bcphilic - endtype - type (tracer_index_type), save :: dmfcs_indx - -! - grid components to be included in the aeropt calculations -!> number of aerosol grid components - integer, save :: num_gridcomp = 0 -!> aerosol grid components - character, allocatable , save :: gridcomp(:)*2 - -!> default full-package setting - integer, parameter :: max_num_gridcomp = 5 -!> data max_gridcomp /'DU', 'BC', 'OC', 'SU', 'SS'/ - character*2 :: max_gridcomp(max_num_gridcomp) - data max_gridcomp /'DU', 'BC', 'OC', 'SU', 'SS'/ - -! GOCART code modification end here (Sarah Lu) -! ------------------------! ! ======================================================================= - +! --------------------------------------------------------------------- ! +! section-5 : module variables for aod diagnostic ! +! --------------------------------------------------------------------- ! !! --- the following are for diagnostic purpose to output aerosol optical depth ! aod from 10 components are grouped into 5 major different species: ! 1:dust (inso,minm,miam,micm,mitr); 2:black carbon (soot) @@ -688,7 +488,6 @@ module module_radiation_aerosols ! public aer_init, aer_update, setaer - ! ================= contains ! ================= @@ -739,7 +538,7 @@ subroutine aer_init & ! ! ! usage: call aer_init ! ! ! -! subprograms called: clim_aerinit, gcrt_aerinit, ! +! subprograms called: clim_aerinit, gocart_aerinit, ! ! wrt_aerlog, set_volcaer, set_spectrum, ! ! ! ! ================================================================== ! @@ -834,14 +633,13 @@ subroutine aer_init & ! --- outputs: & ) -! elseif ( iaermdl == 1 ) then ! gocart-climatology scheme -! elseif ( iaermdl==1 .or. iaermdl==2 ) then ! gocart-clim/prog scheme - -! call gcrt_climinit - -! elseif ( iaermdl == 2 ) then ! gocart-prognostic scheme + elseif ( iaermdl==1 .or. iaermdl==2 ) then ! gocart clim/prog scheme -! call gcrt_aerinit + call gocart_aerinit & +! --- inputs: + & ( solfwv, eirfwv, me & +! --- outputs: + & ) else if ( me == 0 ) then @@ -1962,7 +1760,11 @@ subroutine aer_update & !> -# Call trop_update() to update monthly tropospheric aerosol data. if ( lalwflg .or. laswflg ) then + + if ( iaermdl == 0 .or. iaermdl==5 ) then ! opac-climatology scheme call trop_update + endif + endif !> -# Call volc_update() to update yearly stratospheric volcanic aerosol data. @@ -2294,7 +2096,7 @@ end subroutine aer_update !> @{ !----------------------------------- subroutine setaer & - & ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat, & ! --- inputs + & ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat, & ! --- inputs & IMAX,NLAY,NLP1, lsswr,lslwr, & & aerosw,aerolw & ! --- outputs &, aerodp & @@ -2312,6 +2114,7 @@ subroutine setaer & ! rhlay - layer mean relative humidity IMAX*NLAY ! ! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX ! ! tracer - aerosol tracer concentration IMAX*NLAY*NTRAC ! +! aerfld - prescribed aerosol mixing rat IMAX*NLAY*NTRCAER! ! xlon - longitude of given points in radiance IMAX ! ! ok for both 0->2pi or -pi->+pi ranges ! ! xlat - latitude of given points in radiance IMAX ! @@ -2362,6 +2165,7 @@ subroutine setaer & real (kind=kind_phys), dimension(:), intent(in) :: xlon, xlat, & & slmsk real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer + real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld logical, intent(in) :: lsswr, lslwr @@ -2419,7 +2223,6 @@ subroutine setaer & enddo enddo - if ( .not. (lsswr .or. lslwr) ) then return endif @@ -2495,8 +2298,6 @@ subroutine setaer & !! subroutine computes sw + lw aerosol optical properties for gocart !! aerosol species (merged from fcst and clim fields). -!SARAH -! if ( iaerflg == 1 ) then ! use opac aerosol climatology if ( iaermdl==0 .or. iaermdl==5 ) then ! use opac aerosol climatology call aer_property & @@ -2509,6 +2310,20 @@ subroutine setaer & & aerosw,aerolw,aerodp & & ) +! + elseif ( iaermdl==1 .or. iaermdl==2) then ! use gocart aerosols + + call aer_property_gocart & +! --- inputs: + & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld, & + & alon,alat,slmsk,laersw,laerlw, & + & IMAX,NLAY,NLP1, & +! --- outputs: + & aerosw,aerolw,aerodp & + & ) + endif ! end if_iaerflg_block + + ! --- check print ! do m = 1, NBDSW ! print *,' *** CHECK AEROSOLS PROPERTIES FOR SW BAND =',m, & @@ -2544,21 +2359,6 @@ subroutine setaer & ! print *,' ASYAER:',aerolw(:,k,m,3) ! enddo ! enddo -! SARAH -! elseif ( iaerflg == 2 ) then ! use gocart aerosol scheme - elseif ( iaermdl == 1 ) then ! use gocart aerosol scheme - - call setgocartaer & - -! --- inputs: - & ( alon,alat,prslk,rhlay,dz,hz,NSWLWBD, & - & prsl,tvly,tracer, & - & IMAX,NLAY,NLP1, ivflip, lsswr,lslwr, & -! --- outputs: - & aerosw,aerolw & - & ) - - endif ! end if_iaerflg_block endif ! end if_laswflg_or_lalwflg_block @@ -3267,11 +3067,9 @@ subroutine aer_property & enddo ! --- for diagnostic output (optional) -! if ( lspcaod ) then - do m = 1, NSPC - aerodp(i,m+1) = spcodp(m) - enddo -! endif + do m = 1, NSPC + aerodp(i,m+1) = spcodp(m) + enddo endif ! end if_larsw_block @@ -3613,1499 +3411,824 @@ end subroutine aer_property !----------------------------------- !> @} -! ======================================================================= -! GOCART code modification starts here (Sarah lu) ---------------------! -!! -!! gocart_init : set_aerspc, rd_gocart_clim, rd_gocart_luts, optavg_grt -!! setgocartaer: aeropt_grt, map_aermr - -!> The initialization program for gocart aerosols -!! - determine weight and index for aerosol composition/luts -!! - read in monthly global distribution of gocart aerosols -!! - read and map the tabulated aerosol optical spectral data onto -!! corresponding SW/LW radiation spectral bands. +!> This subroutine is the gocart aerosol initialization +!! program to set up necessary parameters and working arrays. +!>\param solfwv (NWVTOT), solar flux for each individual wavenumber +!! \f$(w/m^2)\f$ +!!\param eirfwv (NWVTIR), IR flux(273k) for each individual wavenumber +!! \f$(w/m^2)\f$ +!!\param me print message control flag !! -!>\param NWVTOT total num of wave numbers used in sw spectrum -!!\param solfwv (NWVTOT), solar flux for each individual -!! wavenumber (w/m2) -!!\param soltot total solar flux for the spectrual range (w/m2) -!!\param NWVTIR total num of wave numbers used in the ir region -!!\param eirfwv (NWVTIR), ir flux(273k) for each individual -!! wavenumber (w/m2) -!!\param NBDSW num of bands calculated for sw aeros opt prop -!!\param NLWBND num of bands calculated for lw aeros opt prop -!!\param NSWLWBD total num of bands calc for sw+lw aeros opt prop -!!\param imon month of the year -!!\param me print message control flag -!!\param raddt -!!\param fdaer !>\section gel_go_ini General Algorithm !! @{ !----------------------------------- - subroutine gocart_init & - & ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv, & ! --- inputs: - & NBDSW,NLWBND,NSWLWBD,imon,me,raddt,fdaer & ! --- outputs: ( none ) + subroutine gocart_aerinit & + & ( solfwv, eirfwv, me & & ) ! ================================================================== ! ! ! -! subprogram : gocart_init ! -! ! -! this is the initialization program for gocart aerosols ! -! ! -! - determine weight and index for aerosol composition/luts ! -! - read in monthly global distribution of gocart aerosols ! -! - read and map the tabulated aerosol optical spectral data ! -! onto corresponding sw/lw radiation spectral bands. ! +! subprogram : gocart_aerinit ! ! ! -! ==================== defination of variables =================== ! +! gocart_aerinit is the gocart aerosol initialization program ! +! to set up necessary parameters and working arrays. ! ! ! ! inputs: ! -! NWVTOT - total num of wave numbers used in sw spectrum ! ! solfwv(NWVTOT) - solar flux for each individual wavenumber (w/m2)! -! soltot - total solar flux for the spectrual range (w/m2)! -! NWVTIR - total num of wave numbers used in the ir region ! ! eirfwv(NWVTIR) - ir flux(273k) for each individual wavenum (w/m2)! -! NBDSW - num of bands calculated for sw aeros opt prop ! -! NLWBND - num of bands calculated for lw aeros opt prop ! -! NSWLWBD - total num of bands calc for sw+lw aeros opt prop! -! imon - month of the year ! ! me - print message control flag ! ! ! -! outputs: (to the module variables) ! +! outputs: (to module variables) ! ! ! ! module variables: ! -! NBDSW - total number of sw spectral bands ! -! wvnum1,wvnum2 (NSWSTR:NSWEND) ! -! - start/end wavenumbers for each of sw bands ! -! NBDLW - total number of lw spectral bands ! -! wvnlw1,wvnlw2 (NBDLW) ! -! - start/end wavenumbers for each of lw bands ! -! NSWLWBD - total number of sw+lw bands used in this version ! -! extrhi_grt - extinction coef for rh-indep aeros KCM1*NSWLWBD ! -! ssarhi_grt - single-scat-alb for rh-indep aeros KCM1*NSWLWBD ! -! asyrhi_grt - asymmetry factor for rh-indep aeros KCM1*NSWLWBD ! -! extrhd_grt - extinction coef for rh-dep aeros KRHLEV*KCM2*NSWLWBD! -! ssarhd_grt - single-scat-alb for rh-dep aeros KRHLEV*KCM2*NSWLWBD! -! asyrhd_grt - asymmetry factor for rh-dep aerosKRHLEV*KCM2*NSWLWBD! -! ctaer - merging coefficients for fcst/clim fields ! -! get_fcst - option to get fcst aerosol fields ! -! get_clim - option to get clim aerosol fields ! -! dm_indx - index for aer spec to be included in aeropt calculations ! -! dmfcs_indx - index for prognostic aerosol fields ! -! psclmg - geos3/4-gocart pressure IMXG*JMXG*KMXG ! -! dmclmg - geos3-gocart aerosol dry mass IMXG*JMXG*KMXG*NMXG! -! or geos4-gocart aerosol mixing ratio ! +! NWVSOL - num of wvnum regions where solar flux is constant ! +! NWVTOT - total num of wave numbers used in sw spectrum ! +! NWVTIR - total num of wave numbers used in the ir region ! +! NSWBND - total number of sw spectral bands ! +! NLWBND - total number of lw spectral bands ! +! NAERBND - number of bands for climatology aerosol data ! +! KCM1 - number of rh independent aeros species ! +! KCM2 - number of rh dependent aeros species ! ! ! ! usage: call gocart_init ! ! ! -! subprograms called: set_aerspc, rd_gocart_clim, ! -! rd_gocart_luts, optavg_grt ! +! subprograms called: rd_gocart_luts, optavg_gocart ! ! ! ! ================================================================== ! implicit none ! --- inputs: - integer, intent(in) :: NWVTOT,NWVTIR,NBDSW,NLWBND,NSWLWBD,imon,me + real (kind=kind_phys), dimension(:) :: solfwv ! one wvn sol flux + real (kind=kind_phys), dimension(:) :: eirfwv ! one wvn ir flux - real (kind=kind_phys), intent(in) :: raddt, fdaer - - real (kind=kind_phys), intent(in) :: solfwv(:),soltot, eirfwv(:) + integer, intent(in) :: me ! --- output: ( none ) ! --- locals: + real (kind=kind_phys), dimension(kaerbndi,kcm1) :: & + & rhidext0_grt, rhidsca0_grt, rhidssa0_grt, rhidasy0_grt + real (kind=kind_phys), dimension(kaerbndd,krhlev,kcm2):: & + & rhdpext0_grt, rhdpsca0_grt, rhdpssa0_grt, rhdpasy0_grt - real (kind=kind_phys), dimension(NBDSW,KAERBND) :: solwaer - real (kind=kind_phys), dimension(NBDSW) :: solbnd - real (kind=kind_phys), dimension(NLWBND,KAERBND) :: eirwaer - real (kind=kind_phys), dimension(NLWBND) :: eirbnd - real (kind=kind_phys) :: sumsol, sumir, fac, tmp, wvs, wve - - integer, dimension(NBDSW) :: nv1, nv2 - integer, dimension(NLWBND) :: nr1, nr2 - - integer :: i, mb, ib, ii, iw, iw1, iw2, ik, ibs, ibe - -!===> ... begin here - -!-------------------------------------------------------------------------- -! (1) determine aerosol specification index and merging coefficients -!-------------------------------------------------------------------------- - - if ( .not. lgrtint ) then - -! --- ... already done aerspc initialization, continue + real (kind=kind_phys), dimension(nswbnd,kaerbndd) :: solwaer + real (kind=kind_phys), dimension(nswbnd) :: solbnd + real (kind=kind_phys), dimension(nlwbnd,kaerbndd) :: eirwaer + real (kind=kind_phys), dimension(nlwbnd) :: eirbnd - continue + real (kind=kind_phys), dimension(nswbnd,kaerbndi) :: solwaer_du + real (kind=kind_phys), dimension(nswbnd) :: solbnd_du + real (kind=kind_phys), dimension(nlwbnd,kaerbndi) :: eirwaer_du + real (kind=kind_phys), dimension(nlwbnd) :: eirbnd_du - else - -! --- ... set aerosol specification index and merging coefficients + integer, dimension(nswbnd) :: nv1, nv2, nv1_du, nv2_du + integer, dimension(nlwbnd) :: nr1, nr2, nr1_du, nr2_du - call set_aerspc(raddt,fdaer) -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) + integer, dimension(kaerbndd) :: iendwv + integer, dimension(kaerbndi) :: iendwv_du + real (kind=kind_phys), dimension(kaerbndd) :: wavelength + real (kind=kind_phys), dimension(kaerbndi) :: wavelength_du + real (kind=kind_phys) :: sumsol, sumir, sumsol_du, sumir_du - endif ! end if_lgrtinit_block + integer :: i, j, k, mb, ib, ii, iix, iw, iw1, iw2 ! -!-------------------------------------------------------------------------- -! (2) read gocart climatological data -!-------------------------------------------------------------------------- - -! --- ... read gocart climatological data, if needed - - if ( get_clim ) then +!===> ... begin here +! +! --- ... invoke gocart aerosol initialization - call rd_gocart_clim -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) + if (KCM /= ntrcaer ) then + print *, 'ERROR in # of gocart aer species',KCM + stop 3000 endif -! -!-------------------------------------------------------------------------- -! (3) read and map the tabulated aerosol optical spectral data -! onto corresponding radiation spectral bands -!-------------------------------------------------------------------------- - - if ( .not. lgrtint ) then +! --- ... aloocate and input aerosol optical data -! --- ... already done optical property interpolation, exit + if ( .not. allocated( extrhi_grt ) ) then + allocate ( extrhi_grt ( kcm1,nswlwbd) ) + allocate ( scarhi_grt ( kcm1,nswlwbd) ) + allocate ( ssarhi_grt ( kcm1,nswlwbd) ) + allocate ( asyrhi_grt ( kcm1,nswlwbd) ) + allocate ( extrhd_grt (krhlev,kcm2,nswlwbd) ) + allocate ( scarhd_grt (krhlev,kcm2,nswlwbd) ) + allocate ( ssarhd_grt (krhlev,kcm2,nswlwbd) ) + allocate ( asyrhd_grt (krhlev,kcm2,nswlwbd) ) + endif - return +! --- ... read tabulated GOCART aerosols optical data - else + call rd_gocart_luts +! --- inputs: (in scope variables, module variables) +! --- outputs: (in scope variables) -! --- ... reset lgrtint +! --- ... convert wavelength to wavenumber +! wavelength and wavelength_du are read-in by rd_gocart_luts - lgrtint = .false. + do i = 1, kaerbndd + iendwv(i) = int(10000. / wavelength(i)) + enddo -! --- ... read tabulated aerosol optical input data - call rd_gocart_luts -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) + do i = 1, kaerbndi + iendwv_du(i) = int(10000. / wavelength_du(i)) + enddo ! --- ... compute solar flux weights and interval indices for mapping ! spectral bands between sw radiation and aerosol data + if ( laswflg ) then solbnd (:) = f_zero - solwaer(:,:) = f_zero + solbnd_du (:)= f_zero + do i=1,nswbnd + do j=1,kaerbndd + solwaer(i,j) = f_zero + enddo + do j=1,kaerbndi + solwaer_du(i,j) = f_zero + enddo + enddo - nv_aod = 1 + do ib = 1, nswbnd + mb = ib + nswstr - 1 + ii = 1 + iix = 1 + iw1 = nint(wvnsw1(mb)) + iw2 = nint(wvnsw2(mb)) - ibs = 1 - ibe = 1 - wvs = wvn_sw1(1) - wve = wvn_sw1(1) - do ib = 2, NBDSW - mb = ib + NSWSTR - 1 - if ( wvn_sw2(mb) >= wvn550 .and. wvn550 >= wvn_sw1(mb) ) then + if ( wvnsw2(mb)>=wvn550 .and. wvn550>=wvnsw1(mb) ) then nv_aod = ib ! sw band number covering 550nm wavelenth endif - if ( wvn_sw1(mb) < wvs ) then - wvs = wvn_sw1(mb) - ibs = ib - endif - if ( wvn_sw1(mb) > wve ) then - wve = wvn_sw1(mb) - ibe = ib - endif - enddo - - do ib = 1, NBDSW - mb = ib + NSWSTR - 1 - ii = 1 - iw1 = nint(wvn_sw1(mb)) - iw2 = nint(wvn_sw2(mb)) - - Lab_swdowhile : do while ( iw1 > iendwv_grt(ii) ) - if ( ii == KAERBND ) exit Lab_swdowhile +! -- for rd-dependent + do while ( iw1 > iendwv(ii) ) + if ( ii == kaerbndd ) exit ii = ii + 1 - enddo Lab_swdowhile - - if ( lmap_new ) then - if (ib == ibs) then - sumsol = f_zero - else - sumsol = -0.5 * solfwv(iw1) - endif - if (ib == ibe) then - fac = f_zero - else - fac = -0.5 - endif - solbnd(ib) = sumsol - else - sumsol = f_zero - endif + enddo + sumsol = f_zero nv1(ib) = ii +! -- for rd-independent + do while ( iw1 > iendwv_du(iix) ) + if ( iix == kaerbndi ) exit + iix = iix + 1 + enddo + sumsol_du = f_zero + nv1_du(ib) = iix + do iw = iw1, iw2 +! -- for rd-dependent solbnd(ib) = solbnd(ib) + solfwv(iw) sumsol = sumsol + solfwv(iw) - if ( iw == iendwv_grt(ii) ) then + if ( iw == iendwv(ii) ) then solwaer(ib,ii) = sumsol - - if ( ii < KAERBND ) then + if ( ii < kaerbndd ) then sumsol = f_zero ii = ii + 1 endif endif + +! -- for rd-independent + solbnd_du(ib) = solbnd_du(ib) + solfwv(iw) + sumsol_du = sumsol_du + solfwv(iw) + + if ( iw == iendwv_du(iix) ) then + solwaer_du(ib,iix) = sumsol_du + if ( iix < kaerbndi ) then + sumsol_du = f_zero + iix = iix + 1 + endif + endif enddo - if ( iw2 /= iendwv_grt(ii) ) then + if ( iw2 /= iendwv(ii) ) then solwaer(ib,ii) = sumsol endif - - if ( lmap_new ) then - tmp = fac * solfwv(iw2) - solwaer(ib,ii) = solwaer(ib,ii) + tmp - solbnd(ib) = solbnd(ib) + tmp + if ( iw2 /= iendwv_du(iix) ) then + solwaer_du(ib,iix) = sumsol_du endif nv2(ib) = ii - - if((me==0) .and. lckprnt) print *,'RAD-nv1,nv2:', & - & ib,iw1,iw2,nv1(ib),iendwv_grt(nv1(ib)), & - & nv2(ib),iendwv_grt(nv2(ib)), & - & 10000./iw1, 10000./iw2 + nv2_du(ib) = iix enddo ! end do_ib_block for sw + endif ! end if_laswflg_block -! --- check the spectral range for the nv_550 band - if((me==0) .and. lckprnt) then - mb = nv_aod + NSWSTR - 1 - iw1 = nint(wvn_sw1(mb)) - iw2 = nint(wvn_sw2(mb)) - print *,'RAD-nv_aod:', & - & nv_aod, iw1, iw2, 10000./iw1, 10000./iw2 - endif -! -! --- ... compute ir flux weights and interval indices for mapping +! --- ... compute lw flux weights and interval indices for mapping ! spectral bands between lw radiation and aerosol data - eirbnd (:) = f_zero - eirwaer(:,:) = f_zero - - ibs = 1 - ibe = 1 - if (NLWBND > 1 ) then - wvs = wvn_lw1(1) - wve = wvn_lw1(1) - do ib = 2, NLWBND - if ( wvn_lw1(ib) < wvs ) then - wvs = wvn_lw1(ib) - ibs = ib - endif - if ( wvn_lw1(ib) > wve ) then - wve = wvn_lw1(ib) - ibe = ib - endif + if ( lalwflg ) then + eirbnd (:) = f_zero + eirbnd_du (:) = f_zero + do i=1,nlwbnd + do j=1,kaerbndd + eirwaer(i,j) = f_zero enddo - endif + do j=1,kaerbndi + eirwaer_du(i,j) = f_zero + enddo + enddo - do ib = 1, NLWBND + do ib = 1, nlwbnd ii = 1 - if ( NLWBND == 1 ) then + iix = 1 + if ( nlwbnd == 1 ) then iw1 = 400 ! corresponding 25 mu iw2 = 2500 ! corresponding 4 mu else - iw1 = nint(wvn_lw1(ib)) - iw2 = nint(wvn_lw2(ib)) + mb = ib + nlwstr - 1 + iw1 = nint(wvnlw1(mb)) + iw2 = nint(wvnlw2(mb)) endif - Lab_lwdowhile : do while ( iw1 > iendwv_grt(ii) ) - if ( ii == KAERBND ) exit Lab_lwdowhile +! -- for rd-dependent + do while ( iw1 > iendwv(ii) ) + if ( ii == kaerbndd ) exit ii = ii + 1 - enddo Lab_lwdowhile - - if ( lmap_new ) then - if (ib == ibs) then - sumir = f_zero - else - sumir = -0.5 * eirfwv(iw1) - endif - if (ib == ibe) then - fac = f_zero - else - fac = -0.5 - endif - eirbnd(ib) = sumir - else - sumir = f_zero - endif + enddo + sumir = f_zero nr1(ib) = ii +! -- for rd-independent + do while ( iw1 > iendwv_du(iix) ) + if ( iix == kaerbndi ) exit + iix = iix + 1 + enddo + sumir_du = f_zero + nr1_du(ib) = iix + do iw = iw1, iw2 +! -- for rd-dependent eirbnd(ib) = eirbnd(ib) + eirfwv(iw) sumir = sumir + eirfwv(iw) - if ( iw == iendwv_grt(ii) ) then + if ( iw == iendwv(ii) ) then eirwaer(ib,ii) = sumir - if ( ii < KAERBND ) then + if ( ii < kaerbndd ) then sumir = f_zero ii = ii + 1 endif endif + +! -- for rd-independent + eirbnd_du(ib) = eirbnd_du(ib) + eirfwv(iw) + sumir_du = sumir_du + eirfwv(iw) + + if ( iw == iendwv_du(iix) ) then + eirwaer_du(ib,iix) = sumir_du + + if ( iix < kaerbndi ) then + sumir_du = f_zero + iix = iix + 1 + endif + endif enddo - if ( iw2 /= iendwv_grt(ii) ) then + if ( iw2 /= iendwv(ii) ) then eirwaer(ib,ii) = sumir endif - - nr2(ib) = ii - - if ( lmap_new ) then - tmp = fac * eirfwv(iw2) - eirwaer(ib,ii) = eirwaer(ib,ii) + tmp - eirbnd(ib) = eirbnd(ib) + tmp + if ( iw2 /= iendwv_du(iix) ) then + eirwaer_du(ib,iix) = sumir_du endif - if(me==0 .and. lckprnt) print *,'RAD-nr1,nr2:', & - & ib,iw1,iw2,nr1(ib),iendwv_grt(nr1(ib)), & - & nr2(ib),iendwv_grt(nr2(ib)), & - & 10000./iw1, 10000./iw2 + nr2(ib) = ii + nr2_du(ib) = iix enddo ! end do_ib_block for lw + endif ! end if_lalwflg_block ! --- compute spectral band mean properties for each species - call optavg_grt -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) - - if(me==0 .and. lckprnt) then - print *, 'RAD -After optavg_grt, sw band info' - do ib = 1, NBDSW - mb = ib + NSWSTR - 1 - print *,'RAD -wvnsw1,wvnsw2: ',ib,wvn_sw1(mb),wvn_sw2(mb) - print *,'RAD -lamda1,lamda2: ',ib,10000./wvn_sw1(mb), & - & 10000./wvn_sw2(mb) - print *,'RAD -extrhi_grt:', extrhi_grt(:,ib) -! do i = 1, KRHLEV - do i = 1, KRHLEV, 10 - print *, 'RAD -extrhd_grt:',i,rhlev_grt(i), & - & extrhd_grt(i,:,ib) - enddo - enddo - print *, 'RAD -After optavg_grt, lw band info' - do ib = 1, NLWBND - ii = NBDSW + ib - print *,'RAD -wvnlw1,wvnlw2: ',ib,wvn_lw1(ib),wvn_lw2(ib) - print *,'RAD -lamda1,lamda2: ',ib,10000./wvn_lw1(ib), & - & 10000./wvn_lw2(ib) - print *,'RAD -extrhi_grt:', extrhi_grt(:,ii) -! do i = 1, KRHLEV - do i = 1, KRHLEV, 10 - print *, 'RAD -extrhd_grt:',i,rhlev_grt(i), & - & extrhd_grt(i,:,ii) - enddo - enddo - endif + call optavg_gocart +! --- inputs: (in-scope variables, module variables) +! --- outputs: (module variables) -! --- ... dealoocate input data arrays no longer needed - deallocate ( iendwv_grt ) - if ( allocated(rhidext0_grt) ) then - deallocate ( rhidext0_grt ) - deallocate ( rhidssa0_grt ) - deallocate ( rhidasy0_grt ) - endif - if ( allocated(rhdpext0_grt) ) then - deallocate ( rhdpext0_grt ) - deallocate ( rhdpssa0_grt ) - deallocate ( rhdpasy0_grt ) - endif - endif ! end if_lgrtinit_block +! --- check print +! if (me == 0) then +! do ib = 1, NSWBND +! mb = ib + NSWSTR - 1 +! print *, ' wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb) +! print *, ' After optavg_gocart, for sw band:',ib +! print *, ' extrhi:', extrhi_grt(:,ib) +! print *, ' scarhi:', scarhi_grt(:,ib) +! print *, ' ssarhi:', ssarhi_grt(:,ib) +! print *, ' asyrhi:', asyrhi_grt(:,ib) +! do i = 1, KRHLEV +! print *, ' extrhd for rhlev:',i +! print *, extrhd_grt(i,:,ib) +! print *, ' scarhd for rhlev:',i +! print *, scarhd_grt(i,:,ib) +! print *, ' ssarhd for rhlev:',i +! print *, ssarhd_grt(i,:,ib) +! print *, ' asyrhd for rhlev:',i +! print *, asyrhd_grt(i,:,ib) +! enddo +! enddo +! print *, ' wvnlw1 :',wvnlw1 +! print *, ' wvnlw2 :',wvnlw2 +! do ib = 1, NLWBND +! ii = NSWBND + ib +! print *,' After optavg_gocart, for lw band:',ib +! print *,' extrhi_grt:', extrhi_grt(:,ii) +! print *,' scarhi_grt:', scarhi_grt(:,ii) +! print *,' ssarhi_grt:', ssarhi_grt(:,ii) +! print *,' asyrhi_grt:', asyrhi_grt(:,ii) +! do i = 1, KRHLEV +! print *,' extrhd for rhlev:',i +! print *, extrhd_grt(i,:,ib) +! print *,' scarhd for rhlev:',i +! print *, scarhd_grt(i,:,ib) +! print *,' ssarhd for rhlev:',i +! print *, ssarhd_grt(i,:,ib) +! print *,' asyrhd for rhlev:',i +! print *, asyrhd_grt(i,:,ib) +! enddo +! enddo +! endif ! ================= contains ! ================= -!> This subroutine determines merging coefficients ctaer; setup aerosol -!! specification. !----------------------------- - subroutine set_aerspc(raddt,fdaer) + subroutine rd_gocart_luts !............................. -! --- inputs: (in scope variables) +! --- inputs: (in scope variables, module variables) ! --- outputs: (in scope variables) ! ==================================================================== ! ! ! -! subprogram: set_aerspc ! -! ! -! determine merging coefficients ctaer; ! -! set up aerosol specification: num_gridcomp, gridcomp, dm_indx, ! -! dmfcs_indx, isoot, iwaso, isuso, issam, isscm ! -! ! -! Aerosol optical properties (ext, ssa, asy) are determined from ! -! NMGX (<=12) aerosol species ! -! ==> DU: dust1 (4 sub-micron bins), dust2, dust3, dust4, dust5 ! -! BC: soot_phobic, soot_philic ! -! OC: waso_phobic, waso_philic ! -! SU: suso (=so4) ! -! SS: ssam (accumulation mode), sscm (coarse mode) ! +! subprogram: rd_gocart_luts ! +! read GMAO pre-tabultaed aerosol optical data for dust, seasalt, ! +! sulfate, black carbon, and organic carbon aerosols ! ! ! -! The current version only supports prognostic aerosols (from GOCART ! -! in-line calculations) and climo aerosols (from GEOS-GOCART runs) ! +! major local variables: ! +! for handling spectral band structures ! +! iendwv - ending wvnum (cm**-1) for each band kaerbndd ! +! iendwv_du - ending wvnum (cm**-1) for each band kaerbndi ! +! for handling optical properties of rh independent species (kcm1) ! +! 1=du001, 2=du002, 3=du003, 4=du004, 5=du005 ! +! rhidext0_grt - extinction coefficient kaerbndi*kcm1 ! +! rhidsca0_grt - scattering coefficient kaerbndi*kcm1 ! +! rhidssa0_grt - single scattering albedo kaerbndi*kcm1 ! +! rhidasy0_grt - asymmetry parameter kaerbndi*kcm1 ! +! for handling optical properties of rh ndependent species (kcm2) ! +! 1=ss001, 2=ss002, 3=ss003, 4=ss004, 5=ss005, 6=so4, ! +! 7=bcphobic, 8=bcphilic, 9=ocphobic, 10=ocphilic ! +! rhdpext0_grt - extinction coefficient kaerbndd*krhlev*kcm2! +! rhdpsca0_grt - scattering coefficient kaerbndd*krhlev*kcm2! +! rhdpssa0_grt - single scattering albedo kaerbndd*krhlev*kcm2! +! rhdpasy0_grt - asymmetry parameter kaerbndd*krhlev*kcm2! +! ! +! usage: call rd_gocart_luts ! ! ! ! ================================================================== ! ! implicit none -! --- inputs: - real (kind=kind_phys), intent(in) :: raddt, fdaer -! --- output: - -! --- local: -! real (kind=kind_phys) :: raddt - integer :: i, indxr - character*2 :: tp, gridcomp_tmp(max_num_gridcomp) - -!! ===> determine ctaer (user specified weight for fcst fields) -! raddt = min(fhswr,fhlwr) / 24. - if( fdaer >= 99999. ) ctaer = f_one - if((fdaer>0.).and.(fdaer<99999.)) ctaer=exp(-raddt/fdaer) - - if(me==0 .and. lckprnt) then - print *, 'RAD -raddt, fdaer,ctaer: ', raddt, fdaer, ctaer - if (ctaer == f_one ) then - print *, 'LU -aerosol fields determined from fcst' - elseif (ctaer == f_zero) then - print *, 'LU -aerosol fields determined from clim' - else - print *, 'LU -aerosol fields determined from fcst/clim' - endif - endif +! --- inputs: (none) +! --- output: (none) -!! ===> determine get_fcst and get_clim -!! if fcst is chosen (ctaer == f_one ), set get_clim to F -!! if clim is chosen (ctaer == f_zero), set get_fcst to F - if ( ctaer == f_one ) get_clim = .false. - if ( ctaer == f_zero ) get_fcst = .false. - -!! ===> determine aerosol species to be included in the calculations -!! of aerosol optical properties (ext, ssa, asy) - -!* If climo option is chosen, the aerosol composition is hardwired -!* to full package. If not, the composition is determined from -!* tracer_config on-the-fly (full package or subset) - lab_if_fcst : if ( get_fcst ) then - -!! use tracer_config to determine num_gridcomp and gridcomp - if ( gfs_phy_tracer%doing_GOCART ) then - if ( gfs_phy_tracer%doing_DU ) then - num_gridcomp = num_gridcomp + 1 - gridcomp_tmp(num_gridcomp) = 'DU' - endif - if ( gfs_phy_tracer%doing_SU ) then - num_gridcomp = num_gridcomp + 1 - gridcomp_tmp(num_gridcomp) = 'SU' - endif - if ( gfs_phy_tracer%doing_SS ) then - num_gridcomp = num_gridcomp + 1 - gridcomp_tmp(num_gridcomp) = 'SS' - endif - if ( gfs_phy_tracer%doing_OC ) then - num_gridcomp = num_gridcomp + 1 - gridcomp_tmp(num_gridcomp) = 'OC' - endif - if ( gfs_phy_tracer%doing_BC ) then - num_gridcomp = num_gridcomp + 1 - gridcomp_tmp(num_gridcomp) = 'BC' - endif +! --- locals: + integer :: iradius, ik, ibeg + integer, parameter :: numspc = 5 ! # of aerosol species + +! - input tabulated aerosol optical spectral data from GSFC + real, dimension(kaerbndd) :: lambda ! wavelength (m) for non-dust + real, dimension(kaerbndi) :: lambda_du ! wavelength (m) for dust + real, dimension(krhlev) :: rh ! relative humidity (fraction) + real, dimension(kaerbndd,krhlev,numspc) :: bext! extinction efficiency (m2/kg) + real, dimension(kaerbndd,krhlev,numspc) :: bsca! scattering efficiency (m2/kg) + real, dimension(kaerbndd,krhlev,numspc) :: g ! asymmetry factor (dimensionless) + real, dimension(kaerbndi,krhlev,numspc) :: bext_du! extinction efficiency (m2/kg) + real, dimension(kaerbndi,krhlev,numspc) :: bsca_du! scattering efficiency (m2/kg) + real, dimension(kaerbndi,krhlev,numspc) :: g_du ! asymmetry factor (dimensionless) ! - if ( num_gridcomp > 0 ) then - allocate ( gridcomp(num_gridcomp) ) - gridcomp(1:num_gridcomp) = gridcomp_tmp(1:num_gridcomp) - else - print *,'ERROR: prognostic aerosols not found,abort',me - stop 1000 - endif - - else ! gfs_phy_tracer%doing_GOCART=F - - print *,'ERROR: prognostic aerosols option off, abort',me - stop 1001 - - endif ! end_if_gfs_phy_tracer%doing_GOCART_if_ - - else lab_if_fcst - -!! set to full package (max_num_gridcomp and max_gridcomp) - num_gridcomp = max_num_gridcomp - allocate ( gridcomp(num_gridcomp) ) - gridcomp(1:num_gridcomp) = max_gridcomp(1:num_gridcomp) - - endif lab_if_fcst - -!! -!! Aerosol specification is determined as such: -!! A. For radiation-aerosol feedback, the specification is based on the aeropt -!! routine from Mian Chin and Hongbin Yu (hydrophobic and hydrophilic for -!! OC/BC; submicron and supermicron for SS, 8-bins (with 4 subgroups for the -!! the submicron bin) for DU, and SO4 for SU) -!! B. For transport, the specification is determined from GOCART in-line module -!! C. For LUTS, (waso, soot, ssam, sscm, suso, dust) is used, based on the -!! the OPAC climo aerosol scheme (implemented by Yu-Tai Hou) - -!!=== determine dm_indx and NMXG - indxr = 0 - dm_indx%waso_phobic = -999 ! OC - dm_indx%soot_phobic = -999 ! BC - dm_indx%ssam = -999 ! SS - dm_indx%suso = -999 ! SU - dm_indx%dust1 = -999 ! DU - do i = 1, num_gridcomp - tp = gridcomp(i) - select case ( tp ) - case ( 'OC') ! consider hydrophobic and hydrophilic - dm_indx%waso_phobic = indxr + 1 - dm_indx%waso_philic = indxr + 2 - indxr = indxr + 2 - case ( 'BC') ! consider hydrophobic and hydrophilic - dm_indx%soot_phobic = indxr + 1 - dm_indx%soot_philic = indxr + 2 - indxr = indxr + 2 - case ( 'SS') ! consider submicron and supermicron - dm_indx%ssam = indxr + 1 - dm_indx%sscm = indxr + 2 - indxr = indxr + 2 - case ( 'SU') ! consider SO4 only - dm_indx%suso = indxr + 1 - indxr = indxr + 1 - case ( 'DU') ! consider all 5 bins - dm_indx%dust1 = indxr + 1 - dm_indx%dust2 = indxr + 2 - dm_indx%dust3 = indxr + 3 - dm_indx%dust4 = indxr + 4 - dm_indx%dust5 = indxr + 5 - indxr = indxr + 5 - case default - print *,'ERROR: aerosol species not supported, abort',me - stop 1002 - end select - enddo -!! - NMXG = indxr ! num of gocart aer spec for opt cal -!! - -!!=== determine dmfcs_indx -!! SS: 5-bins are considered for transport while only two groups -!! (accumulation/coarse modes) are considered for radiation -!! DU: 5-bins are considered for transport while 8 bins (with the -!! submicorn bin exptended to 4 bins) are considered for radiation -!! SU: DMS, SO2, and MSA are not considered for radiation - - if ( get_fcst ) then - if ( gfs_phy_tracer%doing_OC ) then - dmfcs_indx%ocphobic = trcindx ('ocphobic', gfs_phy_tracer) - dmfcs_indx%ocphilic = trcindx ('ocphilic', gfs_phy_tracer) - endif - if ( gfs_phy_tracer%doing_BC ) then - dmfcs_indx%bcphobic = trcindx ('bcphobic', gfs_phy_tracer) - dmfcs_indx%bcphilic = trcindx ('bcphilic', gfs_phy_tracer) - endif - if ( gfs_phy_tracer%doing_SS ) then - dmfcs_indx%ss001 = trcindx ('ss001', gfs_phy_tracer) - dmfcs_indx%ss002 = trcindx ('ss002', gfs_phy_tracer) - dmfcs_indx%ss003 = trcindx ('ss003', gfs_phy_tracer) - dmfcs_indx%ss004 = trcindx ('ss004', gfs_phy_tracer) - dmfcs_indx%ss005 = trcindx ('ss005', gfs_phy_tracer) - endif - if ( gfs_phy_tracer%doing_SU ) then - dmfcs_indx%so4 = trcindx ('so4', gfs_phy_tracer) - endif - if ( gfs_phy_tracer%doing_DU ) then - dmfcs_indx%du001 = trcindx ('du001', gfs_phy_tracer) - dmfcs_indx%du002 = trcindx ('du002', gfs_phy_tracer) - dmfcs_indx%du003 = trcindx ('du003', gfs_phy_tracer) - dmfcs_indx%du004 = trcindx ('du004', gfs_phy_tracer) - dmfcs_indx%du005 = trcindx ('du005', gfs_phy_tracer) - endif - endif + logical :: file_exist + character*50 :: fin, dummy + +! --- read LUTs for dust aerosols + fin='optics_'//gridcomp(1)//'.dat' + inquire (file=trim(fin), exist=file_exist) + if ( file_exist ) then + close(niaercm) + open (unit=niaercm, file=fin, status='OLD') + rewind(niaercm) + else + print *,' Requested luts file ',trim(fin),' not found' + print *,' ** Stopped in rd_gocart_luts ** ' + stop 1220 + endif ! end if_file_exist_block + + iradius = 5 +! read lambda and compute mpwavelength (m) + read(niaercm,'(a40)') dummy + read(niaercm,*) (lambda_du(i), i=1, kaerbndi) +! read rh, relative humidity (fraction) + read(niaercm,'(a40)') dummy + read(niaercm,*) (rh(i), i=1, krhlev) +! read bext (m2 (kg dry mass)-1) + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (bext_du(i,j,k), i=1,kaerbndi) + enddo + enddo +! read bsca (m2 (kg dry mass)-1) + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (bsca_du(i,j,k), i=1, kaerbndi) + enddo + enddo +! read g (dimensionless) + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (g_du(i,j,k), i=1, kaerbndi) + enddo + enddo -!! -!!=== determin KCM, KCM1, KCM2 -!! DU: submicron bin (dust1) contains 4 sub-groups (e.g., hardwire -!! 8 bins for aerosol optical properties luts) -!! OC/BC: while hydrophobic aerosols are rh-independent, the luts -!! for hydrophilic aerosols are used (e.g., use the coeff -!! corresponding to rh=0) -!! - indxr = 1 - isoot = -999 - iwaso = -999 - isuso = -999 - issam = -999 - isscm = -999 - do i = 1, num_gridcomp - tp = gridcomp(i) - if ( tp /= 'DU' ) then !<--- non-dust aerosols - select case ( tp ) - case ( 'OC ') - iwaso = indxr - case ( 'BC ') - isoot = indxr - case ( 'SU ') - isuso = indxr - case ( 'SS ') - issam = indxr - isscm = indxr + 1 - end select - if ( tp /= 'SS' ) then - indxr = indxr + 1 +! fill rhidext0 local arrays for dust aerosols (flip i-index) + do i = 1, kaerbndi ! convert from m to micron + j = kaerbndi -i + 1 ! flip i-index + wavelength_du(j) = 1.e6 * lambda_du(i) + enddo + do k = 1, iradius + do i = 1, kaerbndi + ii = kaerbndi -i + 1 + rhidext0_grt(ii,k) = bext_du(i,1,k) + rhidsca0_grt(ii,k) = bsca_du(i,1,k) + if ( bext_du(i,1,k) /= f_zero) then + rhidssa0_grt(ii,k) = bsca_du(i,1,k)/bext_du(i,1,k) else - indxr = indxr + 2 + rhidssa0_grt(ii,k) = f_one endif - else !<--- dust aerosols - KCM1 = 8 ! num of rh independent aer species - endif - enddo - KCM2 = indxr - 1 ! num of rh dependent aer species - KCM = KCM1 + KCM2 ! total num of aer species - -!! -!! check print starts here - if( me == 0 .and. lckprnt) then - print *, 'RAD -num_gridcomp:', num_gridcomp - print *, 'RAD -gridcomp :', gridcomp(:) - print *, 'RAD -NMXG:', NMXG - print *, 'RAD -dm_indx ===> ' - print *, 'RAD -aerspc: dust1=', dm_indx%dust1 - print *, 'RAD -aerspc: dust2=', dm_indx%dust2 - print *, 'RAD -aerspc: dust3=', dm_indx%dust3 - print *, 'RAD -aerspc: dust4=', dm_indx%dust4 - print *, 'RAD -aerspc: dust5=', dm_indx%dust5 - print *, 'RAD -aerspc: ssam=', dm_indx%ssam - print *, 'RAD -aerspc: sscm=', dm_indx%sscm - print *, 'RAD -aerspc: suso=', dm_indx%suso - print *, 'RAD -aerspc: waso_phobic=',dm_indx%waso_phobic - print *, 'RAD -aerspc: waso_philic=',dm_indx%waso_philic - print *, 'RAD -aerspc: soot_phobic=',dm_indx%soot_phobic - print *, 'RAD -aerspc: soot_philic=',dm_indx%soot_philic - - print *, 'RAD -KCM1 =', KCM1 - print *, 'RAD -KCM2 =', KCM2 - print *, 'RAD -KCM =', KCM - if ( KCM2 > 0 ) then - print *, 'RAD -aerspc: issam=', issam - print *, 'RAD -aerspc: isscm=', isscm - print *, 'RAD -aerspc: isuso=', isuso - print *, 'RAD -aerspc: iwaso=', iwaso - print *, 'RAD -aerspc: isoot=', isoot - endif - - if ( get_fcst ) then - print *, 'RAD -dmfcs_indx ===> ' - print *, 'RAD -trc_du001=',dmfcs_indx%du001 - print *, 'RAD -trc_du002=',dmfcs_indx%du002 - print *, 'RAD -trc_du003=',dmfcs_indx%du003 - print *, 'RAD -trc_du004=',dmfcs_indx%du004 - print *, 'RAD -trc_du005=',dmfcs_indx%du005 - print *, 'RAD -trc_so4 =',dmfcs_indx%so4 - print *, 'RAD -trc_ocphobic=',dmfcs_indx%ocphobic - print *, 'RAD -trc_ocphilic=',dmfcs_indx%ocphilic - print *, 'RAD -trc_bcphobic=',dmfcs_indx%bcphobic - print *, 'RAD -trc_bcphilic=',dmfcs_indx%bcphilic - print *, 'RAD -trc_ss001=',dmfcs_indx%ss001 - print *, 'RAD -trc_ss002=',dmfcs_indx%ss002 - print *, 'RAD -trc_ss003=',dmfcs_indx%ss003 - print *, 'RAD -trc_ss004=',dmfcs_indx%ss004 - print *, 'RAD -trc_ss005=',dmfcs_indx%ss005 - endif - endif -!! check print ends here - - return -! ! - end subroutine set_aerspc - -!----------------------------------- -!> This subroutine reads input gocart aerosol optical data from Mie -!! code calculations. -!----------------------------- - subroutine rd_gocart_luts -!............................. -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) - -! ==================================================================== ! -! subprogram: rd_gocart_luts ! -! read input gocart aerosol optical data from Mie code calculations ! -! ! -! Remarks (Quanhua (Mark) Liu, JCSDA, June 2008) ! -! The LUT is for NCEP selected 61 wave numbers and 6 aerosols ! -! (dust, soot, suso, waso, ssam, and sscm) and 36 aerosol effective ! -! size in microns. ! -! ! -! The LUT is computed using Mie code with a logorithm size ! -! distribution for each of 36 effective sizes. The standard deviation ! -! sigma of the size, and min/max size follows Chin et al. 2000 ! -! For each effective size, it corresponds a relative humidity value. ! -! ! -! The LUT contains the density, sigma, relative humidity, mean mode ! -! radius, effective size, mass extinction coefficient, single ! -! scattering albedo, asymmetry factor, and phase function ! -! ! -! ================================================================== ! -! - implicit none - -! --- inputs: -! --- output: - -! --- locals: - INTEGER, PARAMETER :: NP = 100, NP2 = 2*NP, nWave=100, & - & nAero=6, n_p=36 - INTEGER :: NW, NS, nH, n_bin - real (kind=kind_io8), Dimension( NP2 ) :: Angle, Cos_Angle, & - & Cos_Weight - real (kind=kind_io8), Dimension(n_p,nAero) :: RH, rm, reff - real (kind=kind_io8), Dimension(nWave,n_p,nAero) :: & - & ext0, sca0, asy0 - real (kind=kind_io8), Dimension(NP2,n_p,nWave,nAero) :: ph0 - real (kind=kind_io8) :: wavelength(nWave), density(nAero), & - & sigma(nAero), wave,n_fac,PI,t1,s1,g1 - CHARACTER(len=80) :: AerosolName(nAero) - INTEGER :: i, j, k, l, ij - - character :: aerosol_file*30 - logical :: file_exist - integer :: indx_dust(8) ! map 36 dust bins to gocart size bins - - data aerosol_file /"NCEP_AEROSOL.bin"/ - data AerosolName/ ' Dust ', ' Soot ', ' SUSO ', ' WASO ', & - & ' SSAM ', ' SSCM '/ - -!! 8 dust bins -!! 1 2 3 4 5 6 7 8 -!! .1-.18, .18-.3, .3-.6, 0.6-1.0, 1.0-1.8, 1.8-3, 3-6, 6-10 <-- def -!! 0.1399 0.2399 0.4499 0.8000 1.3994 2.3964 4.4964 7.9887 <-- reff - data indx_dust/4, 8, 12, 18, 21, 24, 30, 36/ - -! PI = acos(-1.d0) - -! -- allocate aerosol optical data - if ( .not. allocated( iendwv_grt ) ) then - allocate ( iendwv_grt (KAERBND) ) - endif - if (.not. allocated(rhidext0_grt) .and. KCM1 > 0 ) then - allocate ( rhidext0_grt(KAERBND,KCM1)) - allocate ( rhidssa0_grt(KAERBND,KCM1)) - allocate ( rhidasy0_grt(KAERBND,KCM1)) - endif - if (.not. allocated(rhdpext0_grt) .and. KCM2 > 0 ) then - allocate ( rhdpext0_grt(KAERBND,KRHLEV,KCM2)) - allocate ( rhdpssa0_grt(KAERBND,KRHLEV,KCM2)) - allocate ( rhdpasy0_grt(KAERBND,KRHLEV,KCM2)) - endif - -! -- read luts - inquire (file = aerosol_file, exist = file_exist) - - if ( file_exist ) then - if(me==0 .and. lckprnt) print *,'RAD -open :',aerosol_file - close (NIAERCM) - open (unit=NIAERCM,file=aerosol_file,status='OLD', & - & action='read',form='UNFORMATTED') - else - print *,' Requested aerosol data file "',aerosol_file, & - & '" not found!', me - print *,' *** Stopped in subroutine RD_GOCART_LUTS !!' - stop 1003 - endif ! end if_file_exist_block - - READ(NIAERCM) (Cos_Angle(i),i=1,NP) - READ(NIAERCM) (Cos_Weight(i),i=1,NP) - READ(NIAERCM) - READ(NIAERCM) - READ(NIAERCM) NW,NS - READ(NIAERCM) - READ(NIAERCM) (wavelength(i),i=1,NW) - -! --- check nAero and NW - if (NW /= KAERBND) then - print *, "Incorrect spectral band, abort ", NW - stop 1004 - endif - -! --- convert wavelength to wavenumber - do i = 1, KAERBND - iendwv_grt(i) = 10000. / wavelength(i) - if(me==0 .and. lckprnt) print *,'RAD -wn,lamda:', & - & i,iendwv_grt(i),wavelength(i) - enddo + rhidasy0_grt(ii,k) = g_du(i,1,k) + enddo + enddo - DO j = 1, nAero - if(me==0 .and. lckprnt) print *,'RAD -read LUTs:', & - & j,AerosolName(j) - READ(NIAERCM) - READ(NIAERCM) - READ(NIAERCM) n_bin, density(j), sigma(j) - READ(NIAERCM) - READ(NIAERCM) (RH(i,j),i=1, n_bin) - READ(NIAERCM) - READ(NIAERCM) (rm(i,j),i=1, n_bin) - READ(NIAERCM) - READ(NIAERCM) (reff(i,j),i=1, n_bin) - -! --- check n_bin - if (n_bin /= KRHLEV ) then - print *, "Incorrect rh levels, abort ", n_bin - stop 1005 - endif +! --- read LUTs for non-dust aerosols + do ib = 2, num_gc ! loop thru SS, SU, BC, OC + fin='optics_'//gridcomp(ib)//'.dat' + inquire (file=trim(fin), exist=file_exist) + if ( file_exist ) then + close(niaercm) + open (unit=niaercm, file=fin, status='OLD') + rewind(niaercm) + else + print *,' Requested luts file ',trim(fin),' not found' + print *,' ** Stopped in rd_gocart_luts ** ' + stop 1222 + endif ! end if_file_exist_block + + ibeg = radius_lower(ib) - kcm1 + iradius = num_radius(ib) + +! read lambda and compute mpwavelength (m) + read(niaercm,'(a40)') dummy + read(niaercm,*) (lambda(i), i=1, kaerbndd) +! read rh, relative humidity (fraction) + read(niaercm,'(a40)') dummy + read(niaercm,*) (rh(i), i=1, krhlev) +! read bext + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (bext(i,j,k), i=1,kaerbndd) + enddo + enddo +! read bsca + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (bsca(i,j,k), i=1, kaerbndd) + enddo + enddo +! read g + do k = 1, iradius + read(niaercm,'(a40)') dummy + do j=1, krhlev + read(niaercm,*) (g(i,j,k), i=1, kaerbndd) + enddo + enddo -! --- read luts - DO k = 1, NW - READ(NIAERCM) wave,(ext0(k,L,j),L=1,n_bin) - READ(NIAERCM) (sca0(k,L,j),L=1,n_bin) - READ(NIAERCM) (asy0(k,L,j),L=1,n_bin) - READ(NIAERCM) (ph0(1:NP2,L,k,j),L=1,n_bin) - END DO - -! --- map luts input to module variables - if (AerosolName(j) == ' Dust ' ) then - if ( KCM1 > 0) then !<-- only if rh independent aerosols are needed - do i = 1, KCM1 - rhidext0_grt(1:KAERBND,i)=ext0(1:KAERBND,indx_dust(i),j) - rhidssa0_grt(1:KAERBND,i)=sca0(1:KAERBND,indx_dust(i),j) - rhidasy0_grt(1:KAERBND,i)=asy0(1:KAERBND,indx_dust(i),j) +! fill rhdpext0 local arrays for non-dust aerosols (flip i-index) + do i = 1, kaerbndd ! convert from m to micron + j = kaerbndd -i + 1 ! flip i-index + wavelength(j) = 1.e6 * lambda(i) + enddo + do k = 1, iradius + ik = ibeg + k - 1 + do i = 1, kaerbndd + ii = kaerbndd -i + 1 + do j = 1, krhlev + rhdpext0_grt(ii,j,ik) = bext(i,j,k) + rhdpsca0_grt(ii,j,ik) = bsca(i,j,k) + if ( bext(i,j,k) /= f_zero) then + rhdpssa0_grt(ii,j,ik) = bsca(i,j,k)/bext(i,j,k) + else + rhdpssa0_grt(ii,j,ik) = f_one + endif + rhdpasy0_grt(ii,j,ik) = g(i,j,k) enddo - endif - else - if ( KCM2 > 0) then !<-- only if rh dependent aerosols are needed - if (AerosolName(j) == ' Soot ') ij = isoot - if (AerosolName(j) == ' SUSO ') ij = isuso - if (AerosolName(j) == ' WASO ') ij = iwaso - if (AerosolName(j) == ' SSAM ') ij = issam - if (AerosolName(j) == ' SSCM ') ij = isscm - if ( ij .ne. -999 ) then - rhdpext0_grt(1:KAERBND,1:KRHLEV,ij) = & - & ext0(1:KAERBND,1:KRHLEV,j) - rhdpssa0_grt(1:KAERBND,1:KRHLEV,ij) = & - & sca0(1:KAERBND,1:KRHLEV,j) - rhdpasy0_grt(1:KAERBND,1:KRHLEV,ij) = & - & asy0(1:KAERBND,1:KRHLEV,j) - endif ! if_ij - endif ! if_KCM2 - endif - END DO + enddo + enddo + + enddo !! ib-loop return !................................... end subroutine rd_gocart_luts !----------------------------------- -! ! -!> This subroutine computes mean aerosols optical properties over each -!! SW/LW radiation spectral band for each of the species components. -!! This program follows GFDL's approach for thick cloud optical property -!! in SW radiation scheme (2000). -!----------------------------- - subroutine optavg_grt -!............................. -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) + +!-------------------------------- + subroutine optavg_gocart +!................................ +! --- inputs: (in-scope variables, module variables) +! --- outputs: (module variables) ! ==================================================================== ! ! ! -! subprogram: optavg_grt ! +! subprogram: optavg_gocart ! ! ! -! compute mean aerosols optical properties over each sw/lw radiation ! +! compute mean aerosol optical properties over each sw radiation ! ! spectral band for each of the species components. This program ! -! follows gfdl's approach for thick cloud opertical property in ! -! sw radiation scheme (2000). ! +! follows optavg routine (in turn follows gfdl's approach for thick ! +! cloud opertical property in sw radiation scheme (2000). ! ! ! ! ==================== defination of variables =================== ! ! ! -! input arguments: ! -! nv1,nv2 (NBDSW) - start/end spectral band indices of aerosol data ! +! major input variables: ! +! nv1,nv2 (nswbnd) - start/end spectral band indices of aerosol data ! ! for each sw radiation spectral band ! -! nr1,nr2 (NLWBND) - start/end spectral band indices of aerosol data ! +! nr1,nr2 (nlwbnd) - start/end spectral band indices of aerosol data ! ! for each ir radiation spectral band ! -! solwaer (NBDSW,KAERBND) ! +! nv1_du,nv2_du(nswbnd) - start/end spectral band indices of aer data! +! for each sw radiation spectral band ! +! nr1_du,nr2_du(nlwbnd) - start/end spectral band indices of aer data! +! for each ir radiation spectral band ! +! solwaer (nswbnd,kaerbndd) ! ! - solar flux weight over each sw radiation band ! ! vs each aerosol data spectral band ! -! eirwaer (NLWBND,KAERBND) ! +! eirwaer (nlwbnd,kaerbndd) ! ! - ir flux weight over each lw radiation band ! ! vs each aerosol data spectral band ! -! solbnd (NBDSW) - solar flux weight over each sw radiation band ! -! eirbnd (NLWBND) - ir flux weight over each lw radiation band ! -! NBDSW - total number of sw spectral bands ! -! NLWBND - total number of lw spectral bands ! -! NSWLWBD - total number of sw+lw spectral bands ! +! solwaer_du (nswbnd,kaerbndi) ! +! - solar flux weight over each sw radiation band ! +! vs each aerosol data spectral band ! +! eirwaer_du (nlwbnd,kaerbndi) ! +! - ir flux weight over each lw radiation band ! +! vs each aerosol data spectral band ! +! solbnd (nswbnd) - solar flux weight over each sw radiation band ! +! eirbnd (nlwbnd) - ir flux weight over each lw radiation band ! +! solbnd_du(nswbnd) - solar flux weight over each sw radiation band ! +! eirbnd_du(nlwbnd) - ir flux weight over each lw radiation band ! +! nswbnd - total number of sw spectral bands ! +! nlwbnd - total number of lw spectral bands ! ! ! -! output arguments: (to module variables) ! +! external module variables: (in physparam) ! +! laswflg - control flag for sw spectral region ! +! lalwflg - control flag for lw spectral region ! +! ! +! output variables: (to module variables) ! ! ! ! ================================================================== ! -! - implicit none ! --- inputs: ! --- output: ! --- locals: - real (kind=kind_phys) :: sumk, sumok, sumokg, sumreft, & + real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft, & & sp, refb, reft, rsolbd, rirbd integer :: ib, nb, ni, nh, nc ! !===> ... begin here - -! --- ... allocate aerosol optical data - if (.not. allocated(extrhd_grt) .and. KCM2 > 0 ) then - allocate ( extrhd_grt(KRHLEV,KCM2,NSWLWBD) ) - allocate ( ssarhd_grt(KRHLEV,KCM2,NSWLWBD) ) - allocate ( asyrhd_grt(KRHLEV,KCM2,NSWLWBD) ) - endif - if (.not. allocated(extrhi_grt) .and. KCM1 > 0 ) then - allocate ( extrhi_grt(KCM1,NSWLWBD) ) - allocate ( ssarhi_grt(KCM1,NSWLWBD) ) - allocate ( asyrhi_grt(KCM1,NSWLWBD) ) - endif ! ! --- ... loop for each sw radiation spectral band - - do nb = 1, NBDSW - rsolbd = f_one / solbnd(nb) - -! --- for rh independent aerosol species - - lab_rhi: if (KCM1 > 0 ) then - do nc = 1, KCM1 - sumk = f_zero - sumok = f_zero - sumokg = f_zero - sumreft = f_zero - - do ni = nv1(nb), nv2(nb) - sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) & - & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) ) - reft = (f_one - sp) / (f_one + sp) - sumreft = sumreft + reft*solwaer(nb,ni) - - sumk = sumk + rhidext0_grt(ni,nc)*solwaer(nb,ni) - sumok = sumok + rhidssa0_grt(ni,nc)*solwaer(nb,ni) & - & * rhidext0_grt(ni,nc) - sumokg = sumokg + rhidssa0_grt(ni,nc)*solwaer(nb,ni) & - & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc) - enddo - - refb = sumreft * rsolbd - - extrhi_grt(nc,nb) = sumk * rsolbd - asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10) - ssarhi_grt(nc,nb) = 4.0*refb & - & / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 ) - - enddo ! end do_nc_block for rh-ind aeros - endif lab_rhi - -! --- for rh dependent aerosols species - - lab_rhd: if (KCM2 > 0 ) then - do nc = 1, KCM2 - do nh = 1, KRHLEV + + if ( laswflg ) then + do nb = 1, nswbnd + rsolbd = f_one / solbnd_du(nb) + do nc = 1, kcm1 ! --- for rh independent aerosol species sumk = f_zero + sums = f_zero sumok = f_zero sumokg = f_zero sumreft = f_zero - do ni = nv1(nb), nv2(nb) - sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) & - & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc))) + do ni = nv1_du(nb), nv2_du(nb) + sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) & + & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) ) reft = (f_one - sp) / (f_one + sp) - sumreft = sumreft + reft*solwaer(nb,ni) - - sumk = sumk + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni) - sumok = sumok + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) & - & * rhdpext0_grt(ni,nh,nc) - sumokg = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) & - & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc) + sumreft = sumreft + reft*solwaer_du(nb,ni) + + sumk = sumk + rhidext0_grt(ni,nc)*solwaer_du(nb,ni) + sums = sums + rhidsca0_grt(ni,nc)*solwaer_du(nb,ni) + sumok = sumok + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) & + & * rhidext0_grt(ni,nc) + sumokg = sumokg + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) & + & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc) enddo refb = sumreft * rsolbd - extrhd_grt(nh,nc,nb) = sumk * rsolbd - asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10) - ssarhd_grt(nh,nc,nb) = 4.0*refb & - & /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2) - enddo ! end do_nh_block - enddo ! end do_nc_block for rh-dep aeros - endif lab_rhd + extrhi_grt(nc,nb) = sumk * rsolbd + scarhi_grt(nc,nb) = sums * rsolbd + asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10) + ssarhi_grt(nc,nb) = 4.0*refb & + & / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 ) + enddo ! end do_nc_block for rh-ind aeros - enddo ! end do_nb_block for sw + rsolbd = f_one / solbnd(nb) + do nc = 1, kcm2 ! --- for rh dependent aerosol species + do nh = 1, krhlev + sumk = f_zero + sums = f_zero + sumok = f_zero + sumokg = f_zero + sumreft = f_zero -! --- ... loop for each lw radiation spectral band + do ni = nv1(nb), nv2(nb) + sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) & + & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc))) + reft = (f_one - sp) / (f_one + sp) + sumreft = sumreft + reft*solwaer(nb,ni) - do nb = 1, NLWBND + sumk = sumk + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni) + sums = sums + rhdpsca0_grt(ni,nh,nc)*solwaer(nb,ni) + sumok = sumok + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) & + & * rhdpext0_grt(ni,nh,nc) + sumokg = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni)& + & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc) + enddo - ib = NBDSW + nb - rirbd = f_one / eirbnd(nb) + refb = sumreft * rsolbd -! --- for rh independent aerosol species + extrhd_grt(nh,nc,nb) = sumk * rsolbd + scarhd_grt(nh,nc,nb) = sums * rsolbd + asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10) + ssarhd_grt(nh,nc,nb) = 4.0*refb & + & /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2) - lab_rhi_lw: if (KCM1 > 0 ) then - do nc = 1, KCM1 - sumk = f_zero - sumok = f_zero - sumokg = f_zero - sumreft = f_zero + enddo ! end do_nh_block + enddo ! end do_nc_block for rh-dep aeros - do ni = nr1(nb), nr2(nb) - sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) & - & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) ) - reft = (f_one - sp) / (f_one + sp) - sumreft = sumreft + reft*eirwaer(nb,ni) - - sumk = sumk + rhidext0_grt(ni,nc)*eirwaer(nb,ni) - sumok = sumok + rhidssa0_grt(ni,nc)*eirwaer(nb,ni) & - & * rhidext0_grt(ni,nc) - sumokg = sumokg + rhidssa0_grt(ni,nc)*eirwaer(nb,ni) & - & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc) - enddo + enddo ! end do_nb_block for sw + endif ! end if_laswflg_block + +! --- ... loop for each lw radiation spectral band - refb = sumreft * rirbd + if ( lalwflg ) then - extrhi_grt(nc,ib) = sumk * rirbd - asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10) - ssarhi_grt(nc,ib) = 4.0*refb & - & / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 ) - enddo ! end do_nc_block for rh-ind aeros - endif lab_rhi_lw + do nb = 1, nlwbnd -! --- for rh dependent aerosols species + ib = nswbnd + nb - lab_rhd_lw: if (KCM2 > 0 ) then - do nc = 1, KCM2 - do nh = 1, KRHLEV + rirbd = f_one / eirbnd_du(nb) + do nc = 1, kcm1 ! --- for rh independent aerosol species sumk = f_zero + sums = f_zero sumok = f_zero sumokg = f_zero sumreft = f_zero - do ni = nr1(nb), nr2(nb) - sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) & - & /(f_one - rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)) ) + do ni = nr1_du(nb), nr2_du(nb) + sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) & + & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) ) reft = (f_one - sp) / (f_one + sp) - sumreft = sumreft + reft*eirwaer(nb,ni) - - sumk = sumk + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni) - sumok = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) & - & * rhdpext0_grt(ni,nh,nc) - sumokg = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) & - & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc) + sumreft = sumreft + reft*eirwaer_du(nb,ni) + + sumk = sumk + rhidext0_grt(ni,nc)*eirwaer_du(nb,ni) + sums = sums + rhidsca0_grt(ni,nc)*eirwaer_du(nb,ni) + sumok = sumok + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) & + & * rhidext0_grt(ni,nc) + sumokg = sumokg + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) & + & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc) enddo refb = sumreft * rirbd - extrhd_grt(nh,nc,ib) = sumk * rirbd - asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10) - ssarhd_grt(nh,nc,ib) = 4.0*refb & - & /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2 ) - enddo ! end do_nh_block - enddo ! end do_nc_block for rh-dep aeros - endif lab_rhd_lw - - enddo ! end do_nb_block for lw - -! - return -!................................ - end subroutine optavg_grt -!-------------------------------- -! -!> This subroutine: -!! - 1. read in aerosol dry mass and surface pressure from GEOS3-GOCART -!! C3.1 2000 monthly dataset or aerosol mixing ratio and surface -!! pressure from GEOS4-GOCART 2000-2007 averaged monthly data set. -!! - 2. compute goes lat/lon array (for horizontal mapping) -!----------------------------------- - subroutine rd_gocart_clim -!................................... -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) - -! ================================================================== ! -! ! -! subprogram: rd_gocart_clim ! -! ! -! 1. read in aerosol dry mass and surface pressure from GEOS3-GOCART ! -! C3.1 2000 monthly data set ! -! or aerosol mixing ratio and surface pressure from GEOS4-GOCART ! -! 2000-2007 averaged monthly data set ! -! 2. compute goes lat/lon array (for horizontal mapping) ! -! ! -! ==================== defination of variables =================== ! -! ! -! inputs arguments: ! -! imon - month of the year ! -! me - print message control flag ! -! ! -! outputs arguments: (to the module variables) ! -! psclmg - pressure (sfc to toa) cb IMXG*JMXG*KMXG ! -! dmclmg - aerosol dry mass/mixing ratio IMXG*JMXG*KMXG*NMXG ! -! geos_rlon - goes longitude deg IMXG ! -! geos_rlat - goes latitude deg JMXG ! -! ! -! usage: call rd_gocart_clim ! -! ! -! program history: ! -! 05/18/2010 --- Lu Add the option to read GEOS4-GOCART climo ! -! ================================================================== ! -! - implicit none - -! --- inputs: -! --- output: + extrhi_grt(nc,ib) = sumk * rirbd + scarhi_grt(nc,ib) = sums * rirbd + asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10) + ssarhi_grt(nc,ib) = 4.0*refb & + & / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 ) -! --- locals: - integer, parameter :: MAXSPC = 5 - real (kind=kind_io4), parameter :: PINT = 0.01 - real (kind=kind_io4), parameter :: EPSQ = 0.0 - - integer :: i, j, k, numspci, ii - integer :: icmp, nrecl, nt1, nt2, nn(MAXSPC) - character :: ymd*6, yr*4, mn*2, tp*2, & - & fname*30, fin*30, aerosol_file*40 - logical :: file_exist - - real (kind=kind_io4), dimension(KMXG) :: sig - real (kind=kind_io4), dimension(IMXG,JMXG) :: ps - real (kind=kind_io4), dimension(IMXG,JMXG,KMXG) :: temp - real (kind=kind_io4), dimension(IMXG,JMXG,KMXG,MAXSPC):: buff - real (kind=kind_phys) :: pstmp - -! Add the following variables for GEOS4-GOCART - real (kind=kind_io4), dimension(KMXG):: hyam, hybm - real (kind=kind_io4) :: p0 - - data yr /'2000'/ !!<=== use 2000 as the climo proxy - -!* sigma_coordinate for GEOS3-GOCART -!* P(i,j,k) = PINT + SIG(k) * (PS(i,j) - PINT) - data SIG / & - & 9.98547E-01,9.94147E-01,9.86350E-01,9.74300E-01,9.56950E-01, & - & 9.33150E-01,9.01750E-01,8.61500E-01,8.11000E-01,7.50600E-01, & - & 6.82900E-01,6.10850E-01,5.37050E-01,4.63900E-01,3.93650E-01, & - & 3.28275E-01,2.69500E-01,2.18295E-01,1.74820E-01,1.38840E-01, & - & 1.09790E-01,8.66900E-02,6.84150E-02,5.39800E-02,4.25750E-02, & - & 3.35700E-02,2.39900E-02,1.36775E-02,5.01750E-03,5.30000E-04 / - -!* hybrid_sigma_pressure_coordinate for GEOS4-GOCART -!* p(i,j,k) = a(k)*p0 + b(k)*ps(i,j) - data hyam/ & - & 0, 0.0062694, 0.02377049, 0.05011813, 0.08278809, 0.1186361, & - & 0.1540329, 0.1836373, 0.2043698, 0.2167788, 0.221193, & - & 0.217729, 0.2062951, 0.1865887, 0.1615213, 0.1372958, & - & 0.1167039, 0.09920014, 0.08432171, 0.06656809, 0.04765031, & - & 0.03382346, 0.0237648, 0.01435208, 0.00659734, 0.002826232, & - & 0.001118959, 0.0004086494, 0.0001368611, 3.750308e-05/ - - data hybm / & - & 0.992555, 0.9642, 0.90556, 0.816375, 0.703815, 0.576585, & - & 0.44445, 0.324385, 0.226815, 0.149165, 0.089375, & - & 0.045865, 0.017485, 0.00348, 0, 0, 0, 0, 0, & - & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 / - - data p0 /1013.25 / - -!===> ... begin here - -! --- allocate and initialize gocart climatological data - if ( .not. allocated (dmclmg) ) then - allocate ( dmclmg(IMXG,JMXG,KMXG,NMXG) ) - allocate ( psclmg(IMXG,JMXG,KMXG) ) - allocate ( molwgt(NMXG) ) - endif - - dmclmg(:,:,:,:) = f_zero - psclmg(:,:,:) = f_zero - molwgt(:) = f_zero + enddo ! end do_nc_block for rh-ind aeros -! --- allocate and initialize geos lat and lon arrays - if ( .not. allocated ( geos_rlon )) then - allocate (geos_rlon(IMXG)) - allocate (geos_rlat(JMXG)) - endif + rirbd = f_one / eirbnd(nb) + do nc = 1, kcm2 ! --- for rh dependent aerosol species + do nh = 1, krhlev + sumk = f_zero + sums = f_zero + sumok = f_zero + sumokg = f_zero + sumreft = f_zero - geos_rlon(:) = f_zero - geos_rlat(:) = f_zero - -! --- compute geos lat and lon arrays - do i = 1, IMXG - geos_rlon(i) = -180. + (i-1)* dltx - end do - do j = 2, JMXG-1 - geos_rlat(j) = -90. + (j-1)* dlty - end do - geos_rlat(1) = -89.5 - geos_rlat(JMXG) = 89.5 - -! --- determine whether GEOS3 or GEOS4 data set is provided - if ( gocart_climo == 'xxxx' ) then - gocart_climo='0000' -! check geos3-gocart climo - aerosol_file = '200001.PS.avg' - inquire (file = aerosol_file, exist = file_exist) - if ( file_exist ) gocart_climo='ver3' -! check geos4-gocart climo - aerosol_file = 'gocart_climo_2000x2007_ps_01.bin' - inquire (file = aerosol_file, exist = file_exist) - if ( file_exist ) gocart_climo='ver4' - endif -! -! -! --- read ps (sfc pressure) and compute 3d pressure field (psclmg) -! - write(mn,'(i2.2)') imon - ymd = yr//mn - aerosol_file = 'null' - if ( gocart_climo == 'ver3' ) then - aerosol_file = ymd//'.PS.avg' - elseif ( gocart_climo == 'ver4' ) then - aerosol_file = 'gocart_climo_2000x2007_ps_'//mn//'.bin' - endif -! - inquire (file = aerosol_file, exist = file_exist) - lab_if_ps : if ( file_exist ) then - - close(NIAERCM) - if ( gocart_climo == 'ver3' ) then - nrecl = 4 * (IMXG * JMXG) - open(NIAERCM, file=trim(aerosol_file), & - & action='read',access='direct',recl=nrecl) - read(NIAERCM, rec=1) ps - do j = 1, JMXG - do i = 1, IMXG - do k = 1, KMXG - pstmp = pint + sig(k) * (ps(i,j) - pint) - psclmg(i,j,k) = 0.1 * pstmp ! convert mb to cb - enddo - enddo - enddo + do ni = nr1(nb), nr2(nb) + sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) & + & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc))) + reft = (f_one - sp) / (f_one + sp) + sumreft = sumreft + reft*eirwaer(nb,ni) - elseif ( gocart_climo == 'ver4' ) then - open(NIAERCM, file=trim(aerosol_file), & - & action='read',status='old', form='unformatted') - read(NIAERCM) ps(:,:) - do j = 1, JMXG - do i = 1, IMXG - do k = 1, KMXG - pstmp = hyam(k)*p0 + hybm(k)*ps(i,j) - psclmg(i,j,k) = 0.1 * pstmp ! convert mb to cb - enddo - enddo - enddo + sumk = sumk + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni) + sums = sums + rhdpsca0_grt(ni,nh,nc)*eirwaer(nb,ni) + sumok = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) & + & * rhdpext0_grt(ni,nh,nc) + sumokg = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) & + & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc) + enddo - endif ! ---- end if_gocart_climo + refb = sumreft * rirbd - else lab_if_ps + extrhd_grt(nh,nc,ib) = sumk * rirbd + scarhd_grt(nh,nc,ib) = sums * rirbd + asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10) + ssarhd_grt(nh,nc,ib) = 4.0*refb & + & /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2) + enddo ! end do_nh_block + enddo ! end do_nc_block for rh-dep aeros - print *,' *** Requested aerosol data file "', & - & trim(aerosol_file), '" not found!' - print *,' *** Stopped in RD_GOCART_CLIM ! ', me - stop 1006 - endif lab_if_ps -! -! --- read aerosol dry mass (g/m3) or mixing ratios (mol/mol,kg/kg) -! - lab_do_icmp : do icmp = 1, num_gridcomp - - tp = gridcomp(icmp) - -! determine aerosol_file - aerosol_file = 'null' - if ( gocart_climo == 'ver3' ) then - if(tp == 'DU') fname='.DU.STD.tv20.g.avg' - if(tp == 'SS') fname='.SS.STD.tv17.g.avg' - if(tp == 'SU') fname='.SU.STD.tv15.g.avg' - if(tp == 'OC') fname='.CC.STD.tv15.g.avg' - if(tp == 'BC') fname='.CC.STD.tv15.g.avg' - aerosol_file=ymd//trim(fname) - elseif ( gocart_climo == 'ver4' ) then - fin = 'gocart_climo_2000x2007_' - if(tp == 'DU') fname=trim(fin)//'du_' - if(tp == 'SS') fname=trim(fin)//'ss_' - if(tp == 'SU') fname=trim(fin)//'su_' - if(tp == 'OC') fname=trim(fin)//'cc_' - if(tp == 'BC') fname=trim(fin)//'cc_' - aerosol_file=trim(fname)//mn//'.bin' - endif - - numspci = 4 - if(tp == 'DU') numspci = 5 - inquire (file=trim(aerosol_file), exist = file_exist) - lab_if_aer: if ( file_exist ) then + enddo ! end do_nb_block for lw + endif ! end if_lalwflg_block ! - close(NIAERCM) - if ( gocart_climo == 'ver3' ) then - nrecl = 4 * numspci * (IMXG * JMXG * KMXG + 3) - open (NIAERCM, file=trim(aerosol_file), & - & action='read',access='direct', recl=nrecl) - read(NIAERCM,rec=1)(nt1,nt2,nn(i),buff(:,:,:,i),i=1,numspci) - - elseif ( gocart_climo == 'ver4' ) then - open (NIAERCM, file=trim(aerosol_file), & - & action='read',status='old', form='unformatted') - do i = 1, numspci - do k = 1, KMXG - read(NIAERCM) temp(:,:,k) - buff(:,:,k,i) = temp(:,:,k) - enddo - enddo - endif - -!!===> fill dmclmg with working array buff - select case ( tp ) - -! fill in DU from DU: du1, du2, du3, du4, du5 - case ('DU' ) - if ( dm_indx%dust1 /= -999) then - do ii = 1, 5 - dmclmg(:,:,:,dm_indx%dust1+ii-1) = buff(:,:,:,ii) - enddo - else - print *, 'ERROR: invalid DU index, abort! ',me - stop 1007 - endif - -! fill in BC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic - case ('BC' ) - if ( dm_indx%soot_phobic /= -999) then - dmclmg(:,:,:,dm_indx%soot_phobic)=buff(:,:,:,1) - dmclmg(:,:,:,dm_indx%soot_philic)=buff(:,:,:,3) - molwgt(dm_indx%soot_phobic) = 12. - molwgt(dm_indx%soot_philic) = 12. - else - print *, 'ERROR: invalid BC index, abort! ',me - stop 1008 - endif - -! fill in SU from SU: dms, so2, so4, msa - case ('SU' ) - if ( dm_indx%suso /= -999) then - dmclmg(:,:,:,dm_indx%suso) = buff(:,:,:,3) - molwgt(dm_indx%suso) = 96. - else - print *, 'ERROR: invalid SU index, abort! ',me - stop 1009 - endif - -! fill in OC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic - case ('OC' ) - if ( dm_indx%waso_phobic /= -999) then - dmclmg(:,:,:,dm_indx%waso_phobic) = 1.4*buff(:,:,:,2) - dmclmg(:,:,:,dm_indx%waso_philic) = 1.4*buff(:,:,:,4) - molwgt(dm_indx%waso_phobic) = 12. - molwgt(dm_indx%waso_philic) = 12. - else - print *, 'ERROR: invalid OC index, abort! ',me - stop 1010 - endif - -! fill in SS from SS: ss1, ss2, ss3, ss4 - case ('SS' ) - if ( dm_indx%ssam /= -999) then - dmclmg(:,:,:,dm_indx%ssam) = buff(:,:,:,1) - dmclmg(:,:,:,dm_indx%sscm) = buff(:,:,:,2) + & - & buff(:,:,:,3)+buff(:,:,:,4) - else - print *, 'ERROR: invalid SS index, abort! ',me - stop 1011 - endif - - case default - - print *, 'ERROR: invalid aerosol species, abort ',tp - stop 1012 - - end select - - else lab_if_aer - print *,' *** Requested aerosol data file "',aerosol_file, & - & '" not found!' - print *,' *** Stopped in RD_GOCART_CLIM ! ', me - stop 1013 - endif lab_if_aer - - enddo lab_do_icmp - + return return !................................... - end subroutine rd_gocart_clim + end subroutine optavg_gocart !----------------------------------- -! + !................................... - end subroutine gocart_init + end subroutine gocart_aerinit !----------------------------------- !! @} -!> This subroutine computes SW + LW aerosol optical properties for -!! gocart aerosol species (merged from fcst and clim fields). -!! -!>\param alon IMAX, longitude of given points in degree -!!\param alat IMAX, latitude of given points in degree -!!\param prslk (IMAX,NLAY), pressure in cb -!!\param rhlay (IMAX,NLAY), layer mean relative humidity -!!\param dz (IMAX,NLAY), layer thickness in m -!!\param hz (IMAX,NLP1), level high in m -!!\param NSWLWBD total number of sw+ir bands for aeros opt prop -!!\param prsl (IMAX,NLAY), layer mean pressure in mb -!!\param tvly (IMAX,NLAY), layer mean virtual temperature in K -!!\param trcly (IMAX,NLAY,NTRAC), layer mean specific tracer in g/g -!!\param IMAX horizontal dimension of arrays -!!\param NLAY,NLP1 vertical dimensions of arrays -!!\param ivflip control flag for direction of vertical index -!!\n =0: index from toa to surface -!!\n =1: index from surface to toa -!!\param lsswr,lslwr logical flag for sw/lw radiation calls -!!\param aerosw (IMAX,NLAY,NBDSW,NF_AESW), aeros opt properties for SW -!!\n (:,:,:,1): optical depth -!!\n (:,:,:,2): single scattering albedo -!!\n (:,:,:,3): asymmetry parameter -!!\param aerolw (IMAX,NLAY,NBDLW,NF_AELW), aeros opt properties for LW -!!\n (:,:,:,1): optical depth -!!\n (:,:,:,2): single scattering albedo -!!\n (:,:,:,3): asymmetry parameter -!>\section gen_setgo General Algorithm -!!@{ +!> This subroutine compute aerosol optical properties for SW +!! and LW radiations. +!!\param prsi (IMAX,NLP1), pressure at interface in mb +!!\param prsl (IMAX,NLAY), layer mean pressure(not used) +!!\param prslk (IMAX,NLAY), exner function=\f$(p/p0)^{rocp}\f$ (not used) +!!\param tvly (IMAX,NLAY), layer virtual temperature (not used) +!!\param rhlay (IMAX,NLAY), layer mean relative humidity +!!\param dz (IMAX,NLAY), layer thickness in m +!!\param hz (IMAX,NLP1), level high in m +!!\param tracer (IMAX,NLAY,NTRAC), aer tracer concentrations +!!\param aerfld (IMAX,NLAY,NTRCAER), aer tracer concentrations +!!\param alon, alat (IMAX), longitude and latitude of given points in degree +!!\param slmsk (IMAX), sea/land mask (sea:0,land:1,sea-ice:2) +!!\param laersw,laerlw logical flag for sw/lw aerosol calculations +!!\param IMAX horizontal dimension of arrays +!!\param NLAY,NLP1 vertical dimensions of arrays +!!\param NSPC num of species for optional aod output fields +!!\param aerosw (IMAX,NLAY,NBDSW,NF_AESW), aeros opt properties for sw +!!\n (:,:,:,1): optical depth +!!\n (:,:,:,2): single scattering albedo +!!\n (:,:,:,3): asymmetry parameter +!!\param aerolw (IMAX,NLAY,NBDLW,NF_AELW), aeros opt properties for lw +!!\n (:,:,:,1): optical depth +!!\n (:,:,:,2): single scattering albedo +!!\n (:,:,:,3): asymmetry parameter +!!\param aerodp (IMAX,NSPC+1), vertically integrated aer-opt-depth +!!\section gel_go_aer_pro General Algorithm +!! @{ !----------------------------------- - subroutine setgocartaer & - & ( alon,alat,prslk,rhlay,dz,hz,NSWLWBD, & ! --- inputs: - & prsl,tvly,trcly, & - & IMAX,NLAY,NLP1, ivflip, lsswr,lslwr, & - & aerosw,aerolw & ! --- outputs: - & ) + subroutine aer_property_gocart & +!................................... +! --- inputs: + & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld, & + & alon,alat,slmsk, laersw,laerlw, & + & imax,nlay,nlp1, & +! --- outputs: + & aerosw,aerolw,aerodp & + & ) ! ================================================================== ! ! ! -! setgocartaer computes sw + lw aerosol optical properties for gocart ! -! aerosol species (merged from fcst and clim fields) ! +! aer_property_gocart maps prescribed gocart aerosol data set onto ! +! model grids, and compute aerosol optical properties for sw and ! +! lw radiations. ! ! ! ! inputs: ! +! prsi - pressure at interface mb IMAX*NLP1 ! +! prsl - layer mean pressure (not used) IMAX*NLAY ! +! prslk - exner function=(p/p0)**rocp (not used) IMAX*NLAY ! +! tvly - layer virtual temperature (not used) IMAX*NLAY ! +! rhlay - layer mean relative humidity IMAX*NLAY ! +! dz - layer thickness m IMAX*NLAY ! +! hz - level high m IMAX*NLP1 ! +! tracer - aer tracer concentrations (not used) IMAX*NLAY*NTRAC! +! aerfld - prescribed aer tracer mixing ratios IMAX*NLAY*NTRCAER! ! alon, alat IMAX ! ! - longitude and latitude of given points in degree ! -! prslk - pressure cb IMAX*NLAY ! -! rhlay - layer mean relative humidity IMAX*NLAY ! -! dz - layer thickness m IMAX*NLAY ! -! hz - level high m IMAX*NLP1 ! -! NSWLWBD - total number of sw+ir bands for aeros opt prop 1 ! -! prsl - layer mean pressure mb IMAX*NLAY ! -! tvly - layer mean virtual temperature k IMAX*NLAY ! -! trcly - layer mean specific tracer g/g IMAX*NLAY*NTRAC! +! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX ! +! laersw,laerlw 1 ! +! - logical flag for sw/lw aerosol calculations ! ! IMAX - horizontal dimension of arrays 1 ! ! NLAY,NLP1-vertical dimensions of arrays 1 ! -! ivflip - control flag for direction of vertical index 1 ! -! =0: index from toa to surface ! -! =1: index from surface to toa ! -! lsswr,lslwr ! -! - logical flag for sw/lw radiation calls 1 ! ! ! ! outputs: ! ! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW! @@ -5116,569 +4239,287 @@ subroutine setgocartaer & ! (:,:,:,1): optical depth ! ! (:,:,:,2): single scattering albedo ! ! (:,:,:,3): asymmetry parameter ! -! tau_gocart - 550nm aeros opt depth IMAX*NLAY*MAX_NUM_GRIDCOMP! +! aerodp - vertically integrated aer-opt-depth IMAX*NSPC+1 ! ! ! ! module parameters and constants: ! -! NBDSW - total number of sw bands for aeros opt prop 1 ! -! NLWBND - total number of ir bands for aeros opt prop 1 ! +! NSWBND - total number of actual sw spectral bands computed ! +! NLWBND - total number of actual lw spectral bands computed ! +! NSWLWBD - total number of sw+lw bands computed ! ! ! -! module variable: (set by subroutine gocart_init) ! -! dmclmg - aerosols dry mass/mixing ratios IMXG*JMXG*KMXG*NMXG ! -! psclmg - pressure cb IMXG*JMXG*KMXG ! +! external module variables: (in physparam) ! +! ivflip - control flag for direction of vertical index ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! ! ! -! usage: call setgocartaer ! +! module variable: (set by subroutine aer_init) ! ! ! -! subprograms called: map_aermr, aeropt_grt ! +! usage: call aer_property_gocart ! ! ! ! ================================================================== ! -! - implicit none ! --- inputs: - integer, intent(in) :: IMAX,NLAY,NLP1,ivflip,NSWLWBD - logical, intent(in) :: lsswr, lslwr + integer, intent(in) :: IMAX, NLAY, NLP1 + logical, intent(in) :: laersw, laerlw - real (kind=kind_phys), dimension(:,:), intent(in) :: prslk, & - & prsl, rhlay, tvly, dz, hz - real (kind=kind_phys), dimension(:), intent(in) :: alon, alat - real (kind=kind_phys), dimension(:,:,:), intent(in) :: trcly + real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, & + & prslk, tvly, rhlay, dz, hz + real (kind=kind_phys), dimension(:), intent(in) :: alon, alat, & + & slmsk + real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer + real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld ! --- outputs: real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: & & aerosw, aerolw + real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp ! --- locals: - real (kind=kind_phys), dimension(NLAY) :: rh1, dz1 - real (kind=kind_phys), dimension(NLAY,NSWLWBD)::tauae,ssaae,asyae - real (kind=kind_phys), dimension(NLAY,max_num_gridcomp) :: & - & tauae_gocart - - real (kind=kind_phys) :: tmp1, tmp2 + real (kind=kind_phys), dimension(nlay,nswlwbd):: tauae,ssaae,asyae + real (kind=kind_phys), dimension(nspc) :: spcodp - integer :: i, i1, i2, j1, j2, k, m, m1, kp - -! prognostic aerosols on gfs grids - real (kind=kind_phys), dimension(:,:,:),allocatable:: aermr,dmfcs - -! aerosol (dry mass) on gfs grids/levels - real (kind=kind_phys), dimension(:,:), allocatable :: & - & dmanl,dmclm, dmclmx - real (kind=kind_phys), dimension(KMXG) :: pstmp, pkstr - real (kind=kind_phys) :: ptop, psfc, tem, plv, tv, rho - -! --- conversion constants - real (kind=kind_phys), parameter :: hdltx = 0.5 * dltx - real (kind=kind_phys), parameter :: hdlty = 0.5 * dlty - -!===> ... begin here -! - if ( .not. allocated(dmanl) ) then - allocate ( dmclmx(KMXG,NMXG) ) - allocate ( dmanl(NLAY,NMXG) ) - allocate ( dmclm(NLAY,NMXG) ) + real (kind=kind_phys),dimension(nlay,kcm) :: aerms + real (kind=kind_phys),dimension(nlay) :: dz1, rh1 + real (kind=kind_phys) :: plv, tv, rho + integer :: i, m, m1, k - allocate ( aermr(IMAX,NLAY,NMXG) ) - allocate ( dmfcs(IMAX,NLAY,NMXG) ) - endif ! -!> -# Call map_aermr() to map input tracer array (trcly) to local -!! tracer array (aermr). - dmfcs(:,:,:) = f_zero - lab_if_fcst : if ( get_fcst ) then - - call map_aermr -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) - - endif lab_if_fcst +!===> ... begin here ! -!> -# Map geos-gocart climo (dmclmg) to gfs grids (dmclm). - lab_do_IMAX : do i = 1, IMAX + lab_do_IMAXg : do i = 1, IMAX - dmclm(:,:) = f_zero - - lab_if_clim : if ( get_clim ) then -! --- map grid in longitude direction - i2 = 1 - j2 = 1 - tmp1 = alon(i) - if (tmp1 > 180.) tmp1 = tmp1 - 360.0 - lab_do_IMXG : do i1 = 1, IMXG - tmp2 = geos_rlon(i1) - if (tmp2 > 180.) tmp2 = tmp2 - 360.0 - if (abs(tmp1-tmp2) <= hdltx) then - i2 = i1 - exit lab_do_IMXG - endif - enddo lab_do_IMXG - -! --- map grid in latitude direction - lab_do_JMXG : do j1 = 1, JMXG - if (abs(alat(i)-geos_rlat(j1)) <= hdlty) then - j2 = j1 - exit lab_do_JMXG - endif - enddo lab_do_JMXG - -! --- update local arrays pstmp and dmclmx - pstmp(:)= psclmg(i2,j2,:)*1000.0 ! cb to Pa - dmclmx(:,:) = dmclmg(i2,j2,:,:) - -! --- map geos-gocart climo (dmclmx) to gfs level (dmclm) - pkstr(:)=fpkap(pstmp(:)) - psfc = pkstr(1) ! pressure at sfc - ptop = pkstr(KMXG) ! pressure at toa - -! --- map grid in verical direction (follow how ozone is mapped -! in radiation_gases routine) +! --- initialize tauae, ssaae, asyae + do m = 1, NSWLWBD do k = 1, NLAY - kp = k ! from sfc to toa - if(ivflip==0) kp = NLAY - k + 1 ! from toa to sfc - tmp1 = prslk(i,kp) - - do m1 = 1, KMXG - 1 ! from sfc to toa - if(tmp1 > pkstr(m1+1) .and. tmp1 <= pkstr(m1)) then - tmp2 = f_one / (pkstr(m1)-pkstr(m1+1)) - tem = (pkstr(m1) - tmp1) * tmp2 - dmclm(kp,:) = tem * dmclmx(m1+1,:)+ & - & (f_one-tem) * dmclmx(m1,:) - endif - enddo - -!* if(tmp1 > psfc) dmclm(kp,:) = dmclmx(1,:) -!* if(tmp1 < ptop) dmclm(kp,:) = dmclmx(KMXG,:) - + tauae(k,m) = f_zero + ssaae(k,m) = f_one + asyae(k,m) = f_zero enddo - endif lab_if_clim -! -! --- compute fcst/clim merged aerosol loading (dmanl) and the -! radiation optical properties (aerosw, aerolw) -! - do k = 1, NLAY + enddo -! --- map global to local arrays (rh1 and dz1) - rh1(k) = rhlay(i,k) - dz1(k) = dz (i,k) +! --- set floor value for aerms (kg/m3) + do k = 1, NLAY + do m = 1, kcm + aerms(k,m) = 1.e-15 + enddo + enddo -! --- convert from mixing ratio to dry mass (g/m3) - plv = 100. * prsl(i,k) ! convert pressure from mb to Pa - tv = tvly(i,k) ! virtual temp in K - rho = plv / (con_rd * tv) ! air density in kg/m3 - if ( get_fcst ) then - do m = 1, NMXG ! mixing ratio (g/g) - dmfcs(i,k,m) = max(1000.*(rho*aermr(i,k,m)),f_zero) - enddo ! m_do_loop - endif - if ( get_clim .and. (gocart_climo == 'ver4') ) then - do m = 1, NMXG - dmclm(k,m)=1000.*dmclm(k,m)*rho !mixing ratio (g/g) - if ( molwgt(m) /= 0. ) then !mixing ratio (mol/mol) - dmclm(k,m)=dmclm(k,m) * (molwgt(m)/con_amd) - endif - enddo ! m_do_loop - endif + do m = 1, nspc + spcodp(m) = f_zero + enddo + do k = 1, NLAY + rh1(k) = rhlay(i,k) ! + dz1(k) = 1000.*dz (i,k) ! thickness converted from km to m + plv = 100.*prsl(i,k) ! convert pressure from mb to Pa + tv = tvly(i,k) ! virtual temp in K + rho = plv / ( con_rd * tv) ! air density in kg/m3 -! --- determine dmanl from dmclm and dmfcs - do m = 1, NMXG - dmanl(k,m)= ctaer*dmfcs(i,k,m) + & - & ( f_one-ctaer)*dmclm(k,m) + do m = 1, KCM + aerms(k,m) = aerfld(i,k,m)*rho ! dry mass (kg/m3) enddo - enddo +! +! --- calculate sw/lw aerosol optical properties for the +! corresponding frequency bands -!> -# Call aeropt_grt() to alculate sw/lw aerosol optical properties -!! for the corresponding frequency bands. + call aeropt +! --- inputs: (in-scope variables) +! --- outputs: (in-scope variables) - call aeropt_grt -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) + enddo ! end_do_k_loop - if ( lsswr ) then +! ---------------------------------------------------------------------- - if ( laswflg ) then +! --- update aerosw and aerolw arrays + if ( laersw ) then - do m = 1, NBDSW - do k = 1, NLAY - aerosw(i,k,m,1) = tauae(k,m) - aerosw(i,k,m,2) = ssaae(k,m) - aerosw(i,k,m,3) = asyae(k,m) - enddo + do m = 1, NBDSW + do k = 1, NLAY + aerosw(i,k,m,1) = tauae(k,m) + aerosw(i,k,m,2) = ssaae(k,m) + aerosw(i,k,m,3) = asyae(k,m) enddo + enddo - else - - aerosw(:,:,:,:) = f_zero - - endif +! --- update diagnostic aod arrays + do k = 1, NLAY + aerodp(i,1) = aerodp(i,1) + tauae(k,nv_aod) + enddo - endif ! end if_lsswr_block + do m = 1, NSPC + aerodp(i,m+1) = spcodp(m) + enddo - if ( lslwr ) then + endif ! end if_larsw_block - if ( lalwflg ) then + if ( laerlw ) then - if ( NLWBND == 1 ) then - m1 = NBDSW + 1 - do m = 1, NBDLW - do k = 1, NLAY - aerolw(i,k,m,1) = tauae(k,m1) - aerolw(i,k,m,2) = ssaae(k,m1) - aerolw(i,k,m,3) = asyae(k,m1) - enddo - enddo - else - do m = 1, NBDLW - m1 = NBDSW + m - do k = 1, NLAY - aerolw(i,k,m,1) = tauae(k,m1) - aerolw(i,k,m,2) = ssaae(k,m1) - aerolw(i,k,m,3) = asyae(k,m1) - enddo + if ( NLWBND == 1 ) then + m1 = NSWBND + 1 + do m = 1, NBDLW + do k = 1, NLAY + aerolw(i,k,m,1) = tauae(k,m1) + aerolw(i,k,m,2) = ssaae(k,m1) + aerolw(i,k,m,3) = asyae(k,m1) enddo - endif - + enddo else - - aerolw(:,:,:,:) = f_zero - + do m = 1, NBDLW + m1 = NSWBND + m + do k = 1, NLAY + aerolw(i,k,m,1) = tauae(k,m1) + aerolw(i,k,m,2) = ssaae(k,m1) + aerolw(i,k,m,3) = asyae(k,m1) + enddo + enddo endif - endif ! end if_lslwr_block - enddo lab_do_IMAX + endif ! end if_laerlw_block + + enddo lab_do_IMAXg ! ================= contains ! ================= -!>\ingroup setaer -!> This subroutine maps input tracer fields (trcly) to local tracer -!! array (aermr). -!----------------------------- - subroutine map_aermr -!............................. -! --- inputs: (in scope variables) -! --- outputs: (in scope variables) - -! ==================================================================== ! -! ! -! subprogram: map_aermr ! -! ! -! map input tracer fields (trcly) to local tracer array (aermr) ! -! ! -! ==================== defination of variables =================== ! -! ! -! input arguments: ! -! IMAX - horizontal dimension of arrays 1 ! -! NLAY - vertical dimensions of arrays 1 ! -! trcly - layer tracer mass mixing ratio g/g IMAX*NLAY*NTRAC! -! output arguments: (to module variables) ! -! aermr - layer aerosol mass mixing ratio g/g IMAX*NLAY*NMXG ! -! ! -! note: ! -! NTRAC is the number of tracers excluding water vapor ! -! NMXG is the number of prognostic aerosol species ! -! ================================================================== ! -! - implicit none - -! --- inputs: -! --- output: - -! --- local: - integer :: i, indx, ii - character :: tp*2 - -! initialize - aermr(:,:,:) = f_zero - ii = 1 !! <---- trcly does not contain q - -! ==> DU: du1 (submicron bins), du2, du3, du4, du5 - if( gfs_phy_tracer%doing_DU ) then - aermr(:,:,dm_indx%dust1) = trcly(:,:,dmfcs_indx%du001-ii) - aermr(:,:,dm_indx%dust2) = trcly(:,:,dmfcs_indx%du002-ii) - aermr(:,:,dm_indx%dust3) = trcly(:,:,dmfcs_indx%du003-ii) - aermr(:,:,dm_indx%dust4) = trcly(:,:,dmfcs_indx%du004-ii) - aermr(:,:,dm_indx%dust5) = trcly(:,:,dmfcs_indx%du005-ii) - endif - -! ==> OC: oc_phobic, oc_philic - if( gfs_phy_tracer%doing_OC ) then - aermr(:,:,dm_indx%waso_phobic) = & - & trcly(:,:,dmfcs_indx%ocphobic-ii) - aermr(:,:,dm_indx%waso_philic) = & - & trcly(:,:,dmfcs_indx%ocphilic-ii) - endif - -! ==> BC: bc_phobic, bc_philic - if( gfs_phy_tracer%doing_BC ) then - aermr(:,:,dm_indx%soot_phobic) = & - & trcly(:,:,dmfcs_indx%bcphobic-ii) - aermr(:,:,dm_indx%soot_philic) = & - & trcly(:,:,dmfcs_indx%bcphilic-ii) - endif - -! ==> SS: ss1, ss2 (submicron bins), ss3, ss4, ss5 - if( gfs_phy_tracer%doing_SS ) then - aermr(:,:,dm_indx%ssam) = trcly(:,:,dmfcs_indx%ss001-ii) & - & + trcly(:,:,dmfcs_indx%ss002-ii) - aermr(:,:,dm_indx%sscm) = trcly(:,:,dmfcs_indx%ss003-ii) & - & + trcly(:,:,dmfcs_indx%ss004-ii) & - & + trcly(:,:,dmfcs_indx%ss005-ii) - endif - -! ==> SU: so4 - if( gfs_phy_tracer%doing_SU ) then - aermr(:,:,dm_indx%suso) = trcly(:,:,dmfcs_indx%so4-ii) - endif - - return -!................................... - end subroutine map_aermr -!----------------------------------- - +!-------------------------------- + subroutine aeropt +!................................ -!>\ingroup setaer -!! This subroutine computes aerosols optical properties in NSWLWBD -!! SW/LW bands. Aerosol distribution at each grid point is composed -!! from up to NMXG aerosol species (from NUM_GRIDCOMP components). -!----------------------------------- - subroutine aeropt_grt -!................................... ! --- inputs: (in scope variables) ! --- outputs: (in scope variables) ! ================================================================== ! ! ! -! subprogram: aeropt_grt ! -! ! -! compute aerosols optical properties in NSWLWBD sw/lw bands. ! -! Aerosol distribution at each grid point is composed from up to ! -! NMXG aerosol species (from NUM_GRIDCOMP components). ! +! compute aerosols optical properties in NSWLWBD bands for gocart ! +! aerosol species ! ! ! ! input variables: ! -! dmanl - aerosol dry mass g/m3 NLAY*NMXG ! ! rh1 - relative humidity % NLAY ! -! dz1 - layer thickness km NLAY ! +! dz1 - layer thickness m NLAY ! +! aerms - aerosol mass concentration kg/m3 NLAY*KCM ! ! NLAY - vertical dimensions - 1 ! -! ivflip - control flag for direction of vertical index ! -! =0: index from toa to surface ! -! =1: index from surface to toa ! ! ! ! output variables: ! -! tauae - aerosol optical depth - NLAY*NSWLWBD ! -! ssaae - aerosol single scattering albedo - NLAY*NSWLWBD ! -! asyae - aerosol asymmetry parameter - NLAY*NSWLWBD ! +! tauae - optical depth - NLAY*NSWLWBD! +! ssaae - single scattering albedo - NLAY*NSWLWBD! +! asyae - asymmetry parameter - NLAY*NSWLWBD! +! aerodp - vertically integrated aer-opt-depth - IMAX*NSPC+1 ! ! ! ! ================================================================== ! -! - implicit none ! --- inputs: ! --- outputs: ! --- locals: - real (kind=kind_phys) :: aerdm - real (kind=kind_phys) :: ext1, ssa1, asy1, ex00, ss00, as00, & - & ex01, ss01, as01, exint - real (kind=kind_phys) :: tau, ssa, asy, & - & sum_tau, sum_ssa, sum_asy - -! --- subgroups for sub-micron dust -! --- corresponds to 0.1-0.18, 0.18-0.3, 0.3-0.6, 0.6-1.0 micron - - real (kind=kind_phys) :: fd(4) - data fd / 0.01053,0.08421,0.25263,0.65263 / - - character :: tp*2 - integer :: icmp, n, kk, ib, ih2, ih1, ii, ij, ijk real (kind=kind_phys) :: drh0, drh1, rdrh - - real (kind=kind_phys) :: qmin !<--lower bound for opt calc - data qmin / 1.e-20 / - -!===> ... begin here - -! --- initialize (assume no aerosols) - tauae = f_zero - ssaae = f_one - asyae = f_zero - - tauae_gocart = f_zero - -!===> ... loop over vertical layers -! - lab_do_layer : do kk = 1, NLAY + real (kind=kind_phys) :: cm, ext01, sca01, asy01, ssa01 + real (kind=kind_phys) :: ext1, asy1, ssa1, sca1 + real (kind=kind_phys) :: sum_tau,sum_asy,sum_ssa,tau,asy,ssa + integer :: ih1, ih2, nbin, ib, ntrc, ktrc ! --- linear interp coeffs for rh-dep species - ih2 = 1 - do while ( rh1(kk) > rhlev_grt(ih2) ) + do while ( rh1(k) > rhlev_grt(ih2) ) ih2 = ih2 + 1 - if ( ih2 > KRHLEV ) exit + if ( ih2 > krhlev ) exit enddo ih1 = max( 1, ih2-1 ) - ih2 = min( KRHLEV, ih2 ) + ih2 = min( krhlev, ih2 ) drh0 = rhlev_grt(ih2) - rhlev_grt(ih1) - drh1 = rh1(kk) - rhlev_grt(ih1) + drh1 = rh1(k) - rhlev_grt(ih1) if ( ih1 == ih2 ) then - rdrh = f_zero + rdrh = f_zero else - rdrh = drh1 / drh0 + rdrh = drh1 / drh0 endif -! --- loop through sw/lw spectral bands - - lab_do_ib : do ib = 1, NSWLWBD - sum_tau = f_zero - sum_ssa = f_zero - sum_asy = f_zero +! --- compute optical properties for each spectral bands + do ib = 1, nswlwbd + + sum_tau = f_zero + sum_ssa = f_zero + sum_asy = f_zero + +! --- determine tau, ssa, asy for dust aerosols + ext1 = f_zero + asy1 = f_zero + sca1 = f_zero + ssa1 = f_zero + do m = 1, kcm1 + cm = max(aerms(k,m),0.0) * dz1(k) + ext1 = ext1 + cm*extrhi_grt(m,ib) + sca1 = sca1 + cm*scarhi_grt(m,ib) + ssa1 = ssa1 + cm*extrhi_grt(m,ib) * ssarhi_grt(m,ib) + asy1 = asy1 + cm*scarhi_grt(m,ib) * asyrhi_grt(m,ib) + enddo ! m-loop + tau = ext1 + if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1) + if (sca1 > f_zero) asy=min(f_one, asy1/sca1) + +! --- update aod from individual species + if ( ib==nv_aod ) then + spcodp(1) = spcodp(1) + tau + endif +! --- update sum_tau, sum_ssa, sum_asy + sum_tau = sum_tau + tau + sum_ssa = sum_ssa + tau * ssa + sum_asy = sum_asy + tau * ssa * asy -! --- loop through aerosol grid components - lab_do_icmp : do icmp = 1, NUM_GRIDCOMP +! --- determine tau, ssa, asy for non-dust aerosols + do ntrc = 2, nspc ext1 = f_zero - ssa1 = f_zero asy1 = f_zero - - tp = gridcomp(icmp) - - select case ( tp ) - -! -- dust aerosols: no humidification effect - case ( 'DU') - do n = 1, KCM1 - - if (n <= 4) then - aerdm = dmanl(kk,dm_indx%dust1) * fd(n) - else - aerdm = dmanl(kk,dm_indx%dust1+n-4 ) - endif - - if (aerdm < qmin) aerdm = f_zero - ex00 = extrhi_grt(n,ib)*(1000.*dz1(kk))*aerdm - ss00 = ssarhi_grt(n,ib) - as00 = asyrhi_grt(n,ib) - ext1 = ext1 + ex00 - ssa1 = ssa1 + ex00 * ss00 - asy1 = asy1 + ex00 * ss00 * as00 - - enddo - -! -- suso aerosols: with humidification effect - case ( 'SU') - ij = isuso - exint = extrhd_grt(ih1,ij,ib) & - & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib)) - ss00 = ssarhd_grt(ih1,ij,ib) & - & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib)) - as00 = asyrhd_grt(ih1,ij,ib) & - & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib)) - - aerdm = dmanl(kk, dm_indx%suso) - if (aerdm < qmin) aerdm = f_zero - ex00 = exint*(1000.*dz1(kk))*aerdm - ext1 = ex00 - ssa1 = ex00 * ss00 - asy1 = ex00 * ss00 * as00 - -! -- seasalt aerosols: with humidification effect - case ( 'SS') - do n = 1, 2 !<---- ssam, sscm - ij = issam + (n-1) - exint = extrhd_grt(ih1,ij,ib) & - & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib)) - ss00 = ssarhd_grt(ih1,ij,ib) & - & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib)) - as00 = asyrhd_grt(ih1,ij,ib) & - & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib)) - - aerdm = dmanl(kk, dm_indx%ssam+n-1) - if (aerdm < qmin) aerdm = f_zero - ex00 = exint*(1000.*dz1(kk))*aerdm - ext1 = ext1 + ex00 - ssa1 = ssa1 + ex00 * ss00 - asy1 = asy1 + ex00 * ss00 * as00 - - enddo - -! -- organic carbon/black carbon: -! using 'waso' and 'soot' for hydrophilic OC and BC -! using 'waso' and 'soot' at RH=0 for hydrophobic OC and BC - case ( 'OC', 'BC') - if(tp == 'OC') then - ii = dm_indx%waso_phobic - ij = iwaso - else - ii = dm_indx%soot_phobic - ij = isoot - endif - -! --- hydrophobic - aerdm = dmanl(kk, ii) - if (aerdm < qmin) aerdm = f_zero - ex00 = extrhd_grt(1,ij,ib)*(1000.*dz1(kk))*aerdm - ss00 = ssarhd_grt(1,ij,ib) - as00 = asyrhd_grt(1,ij,ib) -! --- hydrophilic - aerdm = dmanl(kk, ii+1) - if (aerdm < qmin) aerdm = f_zero - exint = extrhd_grt(ih1,ij,ib) & - & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib)) - ex01 = exint*(1000.*dz1(kk))*aerdm - ss01 = ssarhd_grt(ih1,ij,ib) & - & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib)) - as01 = asyrhd_grt(ih1,ij,ib) & - & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib)) - - ext1 = ex00 + ex01 - ssa1 = (ex00 * ss00) + (ex01 * ss01) - asy1 = (ex00 * ss00 * as00) + (ex01 * ss01 * as01) - - end select - -! --- determine tau, ssa, asy for each grid component + sca1 = f_zero + ssa1 = f_zero + ktrc = trc_to_aod(ntrc) + do nbin = 1, num_radius(ntrc) + m1 = radius_lower(ntrc) + nbin - 1 + m = m1 - num_radius(1) ! exclude dust aerosols + cm = max(aerms(k,m1),0.0) * dz1(k) + ext01 = extrhd_grt(ih1,m,ib) + & + & rdrh * (extrhd_grt(ih2,m,ib)-extrhd_grt(ih1,m,ib)) + sca01 = scarhd_grt(ih1,m,ib) + & + & rdrh * (scarhd_grt(ih2,m,ib)-scarhd_grt(ih1,m,ib)) + ssa01 = ssarhd_grt(ih1,m,ib) + & + & rdrh * (ssarhd_grt(ih2,m,ib)-ssarhd_grt(ih1,m,ib)) + asy01 = asyrhd_grt(ih1,m,ib) + & + & rdrh * (asyrhd_grt(ih2,m,ib)-asyrhd_grt(ih1,m,ib)) + ext1 = ext1 + cm*ext01 + sca1 = sca1 + cm*sca01 + ssa1 = ssa1 + cm*ext01 * ssa01 + asy1 = asy1 + cm*sca01 * asy01 + enddo ! end_do_nbin_loop tau = ext1 - if (ext1 > f_zero) ssa=min(f_one,ssa1/ext1) - if (ssa1 > f_zero) asy=min(f_one,asy1/ssa1) - -! --- save tau at 550 nm for each grid component - if ( ib == nv_aod ) then - do ijk = 1, max_num_gridcomp - if ( tp == max_gridcomp(ijk) ) then - tauae_gocart(kk,ijk) = tau - endif - enddo + if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1) + if (sca1 > f_zero) asy=min(f_one, asy1/sca1) +! --- update aod from individual species + if ( ib==nv_aod ) then + spcodp(ktrc) = spcodp(ktrc) + tau endif - ! --- update sum_tau, sum_ssa, sum_asy sum_tau = sum_tau + tau sum_ssa = sum_ssa + tau * ssa sum_asy = sum_asy + tau * ssa * asy - - enddo lab_do_icmp - + enddo ! end_do_ntrc_loop ! --- determine total tau, ssa, asy for aerosol mixture - tauae(kk,ib) = sum_tau - if (sum_tau > f_zero) ssaae(kk,ib) = sum_ssa / sum_tau - if (sum_ssa > f_zero) asyae(kk,ib) = sum_asy / sum_ssa - - enddo lab_do_ib - - enddo lab_do_layer + tauae(k,ib) = sum_tau + if (sum_tau > f_zero) ssaae(k,ib) = sum_ssa / sum_tau + if (sum_ssa > f_zero) asyae(k,ib) = sum_asy / sum_ssa + enddo ! end_do_ib_loop ! return -!................................... - end subroutine aeropt_grt -!-------------------------------- - !................................ - end subroutine setgocartaer + end subroutine aeropt !-------------------------------- + +!................................... + end subroutine aer_property_gocart +!----------------------------------- !! @} ! -! GOCART code modification end here (Sarah Lu) ------------------------! ! ======================================================================= !..........................................! From a3b141c75fdb5b4363c8c933917897b1b44ccec1 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 27 Jan 2020 15:11:53 +0000 Subject: [PATCH 09/21] fixed too much high level cloud for iccn==2 --- ccpp/physics | 2 +- gfsphysics/physics/m_micro_driver.F90 | 24 +++++++++++++----------- gfsphysics/physics/micro_mg3_0.F90 | 6 +++--- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/ccpp/physics b/ccpp/physics index 47ecb07ac..11821ddcb 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 47ecb07ac0e7dfee9537fa5a013b7bf1e9ed9b7c +Subproject commit 11821ddcbe34ff4ed53950f14348911914a56c7e diff --git a/gfsphysics/physics/m_micro_driver.F90 b/gfsphysics/physics/m_micro_driver.F90 index 203fe240e..839ef0c4a 100644 --- a/gfsphysics/physics/m_micro_driver.F90 +++ b/gfsphysics/physics/m_micro_driver.F90 @@ -830,17 +830,19 @@ subroutine m_micro_driver(im, ix, lm, flipv, dt_i & ! if(temp(i,k) > T_ICE_ALL) SC_ICE(i,k) = 1.0 ! if(temp(i,k) > TICE) SC_ICE(i,k) = rhc(i,k) ! - if(temp(i,k) < T_ICE_ALL) then -! SC_ICE(i,k) = max(SC_ICE(I,k), 1.2) - SC_ICE(i,k) = max(SC_ICE(I,k), 1.5) - elseif(temp(i,k) > TICE) then - SC_ICE(i,k) = rhc(i,k) - else -! SC_ICE(i,k) = 1.0 -! tx1 = max(SC_ICE(I,k), 1.2) - tx1 = max(SC_ICE(I,k), 1.5) - SC_ICE(i,k) = ((tice-temp(i,k))*tx1 + (temp(i,k)-t_ice_all)*rhc(i,k)) & - * t_ice_denom + if(iccn == 0) then + if(temp(i,k) < T_ICE_ALL) then +! SC_ICE(i,k) = max(SC_ICE(I,k), 1.2) + SC_ICE(i,k) = max(SC_ICE(I,k), 1.5) + elseif(temp(i,k) > TICE) then + SC_ICE(i,k) = rhc(i,k) + else +! SC_ICE(i,k) = 1.0 +! tx1 = max(SC_ICE(I,k), 1.2) + tx1 = max(SC_ICE(I,k), 1.5) + SC_ICE(i,k) = ((tice-temp(i,k))*tx1 + & + (temp(i,k)-t_ice_all)*rhc(i,k))* t_ice_denom + endif endif if (iccn.ne.1) then CDNC_NUC(I,k) = npccninr8(k) diff --git a/gfsphysics/physics/micro_mg3_0.F90 b/gfsphysics/physics/micro_mg3_0.F90 index db0870bee..386a4794a 100644 --- a/gfsphysics/physics/micro_mg3_0.F90 +++ b/gfsphysics/physics/micro_mg3_0.F90 @@ -594,8 +594,7 @@ subroutine micro_mg_tend ( & real(r8), intent(in) :: icecldf(mgncol,nlev) ! ice cloud fraction (no units) real(r8), intent(in) :: qsatfac(mgncol,nlev) ! subgrid cloud water saturation scaling factor (no units) logical, intent(in) :: lprnt - integer, intent(in) :: iccn - + integer, intent(in) :: iccn ! used for scavenging ! Inputs for aerosol activation @@ -1516,6 +1515,8 @@ subroutine micro_mg_tend ( & do i=1,mgncol if (t(i,k) < icenuct) then ncai(i,k) = naai(i,k)*rho(i,k) + ncai(i,k) = min(ncai(i,k), 710.0e3_r8) + naai(i,k) = ncai(i,k)*rhoinv(i,k) else naai(i,k) = zero ncai(i,k) = zero @@ -2834,7 +2835,6 @@ subroutine micro_mg_tend ( & qctend(i,k) = qctend(i,k)+ & (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- & psacws(i,k)-bergs(i,k)-qmultg(i,k)-psacwg(i,k)-pgsacw(i,k))*lcldm(i,k)-berg(i,k) - if (do_cldice) then ! qitend(i,k) = qitend(i,k) + & ! (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+(-prci(i,k)- & From dc2aa92bdea2d07a2bf7ebda756504dd3954f1ef Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Thu, 13 Feb 2020 12:44:30 -0500 Subject: [PATCH 10/21] changed ntrcaer in rad_aero to ntrcaerm --- ccpp/suites/suite_FV3_GFS_v16_csawmg.xml | 2 +- gfsphysics/physics/radiation_aerosols.f | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml index 147d00891..9e2dbed22 100644 --- a/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml +++ b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml @@ -1,6 +1,6 @@ - + diff --git a/gfsphysics/physics/radiation_aerosols.f b/gfsphysics/physics/radiation_aerosols.f index 339b991f0..45a909ca8 100644 --- a/gfsphysics/physics/radiation_aerosols.f +++ b/gfsphysics/physics/radiation_aerosols.f @@ -169,7 +169,7 @@ module module_radiation_aerosols ! use module_radlw_parameters, only : NBDLW, wvnlw1, wvnlw2 ! use funcphys, only : fpkap - use aerclm_def, only : ntrcaer + use aerclm_def, only : ntrcaerm ! implicit none @@ -3499,7 +3499,7 @@ subroutine gocart_aerinit & ! --- ... invoke gocart aerosol initialization - if (KCM /= ntrcaer ) then + if (KCM /= ntrcaerm ) then print *, 'ERROR in # of gocart aer species',KCM stop 3000 endif From 5f82a9502e788c817181c5e28789310ef70dffd6 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Thu, 19 Mar 2020 16:55:00 -0400 Subject: [PATCH 11/21] commited on MG3_v1 on 03/19/2020 --- ccpp/physics | 2 +- gfsphysics/physics/aerinterp.f90 | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 11821ddcb..6466b9105 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 11821ddcbe34ff4ed53950f14348911914a56c7e +Subproject commit 6466b9105bdad7830467c5035c31b16205d7f795 diff --git a/gfsphysics/physics/aerinterp.f90 b/gfsphysics/physics/aerinterp.f90 index 8d35efbcf..6a0eadb99 100644 --- a/gfsphysics/physics/aerinterp.f90 +++ b/gfsphysics/physics/aerinterp.f90 @@ -162,7 +162,13 @@ SUBROUTINE read_aerdata (me, master, iflip, idate ) endif do i = 1, hmx aerin(i+hmx,j,k,ii,imon) = 1.d0*buffx(i,j,klev,1) + if(aerin(i+hmx,j,k,ii,imon)<0.or.aerin(i+hmx,j,k,ii,imon)>1.) then + aerin(i+hmx,j,k,ii,imon) = 0. + end if aerin(i,j,k,ii,imon) = 1.d0*buffx(i+hmx,j,klev,1) + if(aerin(i,j,k,ii,imon)<0.or.aerin(i,j,k,ii,imon)>1.) then + aerin(i,j,k,ii,imon) = 0. + end if enddo !i-loop (lon) enddo !k-loop (lev) enddo !j-loop (lat) From c4bc462c80af136d8dfeb8436c7be57fbbb29ae8 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 23 Mar 2020 13:24:47 -0400 Subject: [PATCH 12/21] changes made in gitmodules --- .gitmodules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index d253f6966..a4bf59883 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,5 +8,5 @@ branch = master [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/NCAR/ccpp-physics + url = https://github.com/AnningCheng-NOAA/ccpp-physics branch = master From 1811200986db1dc6c0ac1f95b2c7673ad433011b Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 24 Mar 2020 16:01:21 -0400 Subject: [PATCH 13/21] update fv3 dycore and ccpp/framework to current develop version --- atmos_cubed_sphere | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere index 0e84f88b4..371a29afb 160000 --- a/atmos_cubed_sphere +++ b/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit 0e84f88b494b9e0a4097da50abe6b143330e8a2f +Subproject commit 371a29afbf813357dd93647cac0cbcd44db2ab20 From f5c9ac3fe7bdc1dcd2e52eccf2fb8faed547adf8 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Tue, 24 Mar 2020 16:46:11 -0400 Subject: [PATCH 14/21] Update atmos_cubed_sphere and ccpp/framework --- ccpp/framework | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/framework b/ccpp/framework index 28fdedec6..d32b965b1 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit 28fdedec637b3a11a90aa3d816aff9dd56a9dc5e +Subproject commit d32b965b11882a42d9db522dc13823b7720b63aa From 0d5a9723ce37f1cc210c1f99a2721152d2b55fbe Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Wed, 25 Mar 2020 16:58:04 -0400 Subject: [PATCH 15/21] update suite_FV3_GFS_v16_csawmg.xml to fix compilation issue such as use machine and py program error caused by dcyc2t3_post --- ccpp/physics | 2 +- ccpp/suites/suite_FV3_GFS_v16_csawmg.xml | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/ccpp/physics b/ccpp/physics index bdc185afa..47713ac7c 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit bdc185afa0467fddfa33e14b7cc0843f29fc9ea7 +Subproject commit 47713ac7cd1cd75cc9f74ca6dc109f8af29a0d5e diff --git a/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml index 9e2dbed22..f1522d7fa 100644 --- a/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml +++ b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml @@ -49,7 +49,6 @@ GFS_surface_composites_post - dcyc2t3_post sfc_diag sfc_diag_post GFS_surface_generic_post @@ -59,6 +58,7 @@ GFS_GWD_generic_pre cires_ugwp cires_ugwp_post + GFS_GWD_generic_post rayleigh_damp GFS_suite_stateout_update ozphys_2015 @@ -72,7 +72,6 @@ GFS_DCNV_generic_post GFS_SCNV_generic_pre samfshalcnv - samfshalcnv_post GFS_SCNV_generic_post GFS_suite_interstitial_4 cnvc90 From 2de35c86da2054bb8fdb0f46c273b19835681987 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Sun, 29 Mar 2020 01:36:56 +0000 Subject: [PATCH 16/21] merge with develop 3/28/2020 --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 47713ac7c..84daf0cb3 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 47713ac7cd1cd75cc9f74ca6dc109f8af29a0d5e +Subproject commit 84daf0cb31446e6998dbd1bc40d6bcf6d86e826a From 4a7650c85eddfc46fdba4c6853fec3da63cef711 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 30 Mar 2020 09:40:46 -0400 Subject: [PATCH 17/21] keep my changes before merge with fv3atm --- ccpp/physics | 2 +- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 93 ++++++++++----------- gfsphysics/GFS_layer/GFS_typedefs.F90 | 46 ++++------ gfsphysics/GFS_layer/GFS_typedefs.meta | 35 -------- gfsphysics/physics/GFS_debug.F90 | 24 +++++- gfsphysics/physics/gcycle.F90 | 12 +-- 6 files changed, 86 insertions(+), 126 deletions(-) diff --git a/ccpp/physics b/ccpp/physics index 84daf0cb3..47713ac7c 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 84daf0cb31446e6998dbd1bc40d6bcf6d86e826a +Subproject commit 47713ac7cd1cd75cc9f74ca6dc109f8af29a0d5e diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 5725c8ba1..1929f030b 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -544,9 +544,6 @@ subroutine GFS_physics_driver & doms, psautco_l, prautco_l, ocalnirbm_cpl, ocalnirdf_cpl, & ocalvisbm_cpl, ocalvisdf_cpl, dtzm, temrain1, t2mmp, q2mp, & psaur_l, praur_l, & -!--- coupling inputs for physics - dtsfc_cice, dqsfc_cice, dusfc_cice, dvsfc_cice, & -! dtsfc_cice, dqsfc_cice, dusfc_cice, dvsfc_cice, ulwsfc_cice, & !--- for CS-convection wcbmax @@ -678,16 +675,15 @@ subroutine GFS_physics_driver & real :: pshltr,QCQ,rh02 real(kind=kind_phys), allocatable, dimension(:,:) :: den - !! Initialize local variables (mainly for debugging purposes, because the - !! corresponding variables Interstitial(nt)%... are reset to zero every time); - !! these variables are only modified over parts of the entire domain (related - !! to land surface mask etc.) + !! Initialize local variables (for debugging purposes only, + !! because the corresponding variables Interstitial(nt)%... + !! are reset to zero every time). !snowmt = 0. !gamq = 0. !gamt = 0. !gflx = 0. !hflx = 0. - ! + !! Strictly speaking, this is not required. But when !! hunting for bit-for-bit differences, doing the same as !! in GFS_suite_stateout_reset makes life a lot easier. @@ -911,10 +907,6 @@ subroutine GFS_physics_driver & ! --- set initial quantities for stochastic physics deltas if (Model%do_sppt) then Tbd%dtdtr = zero - do i=1,im - Tbd%drain_cpl(i) = Coupling%rain_cpl (i) - Tbd%dsnow_cpl(i) = Coupling%snow_cpl (i) - enddo endif ! mg, sfc-perts @@ -1113,14 +1105,6 @@ subroutine GFS_physics_driver & do i=1,im islmsk_cice(i) = nint(Coupling%slimskin_cpl(i)) flag_cice(i) = (islmsk_cice(i) == 4) - - if (flag_cice(i)) then -! ulwsfc_cice(i) = Coupling%ulwsfcin_cpl(i) - dusfc_cice(i) = Coupling%dusfcin_cpl(i) - dvsfc_cice(i) = Coupling%dvsfcin_cpl(i) - dtsfc_cice(i) = Coupling%dtsfcin_cpl(i) - dqsfc_cice(i) = Coupling%dqsfcin_cpl(i) - endif enddo endif !*## CCPP ## @@ -1145,9 +1129,8 @@ subroutine GFS_physics_driver & endif endif if (fice(i) < one) then - wet(i)=.true. !some open ocean/lake water exists + wet(i)=.true. ! some open ocean/lake water exists if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) - end if else fice(i) = zero @@ -1659,13 +1642,14 @@ subroutine GFS_physics_driver & sbsno(i) = zero snowc(i) = zero snohf(i) = zero + !## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run Diag%zlvl(i) = Statein%phil(i,1) * onebg Diag%smcwlt2(i) = zero Diag%smcref2(i) = zero - wind(i) = max(sqrt(Statein%ugrs(i,1)*Statein%ugrs(i,1) + & Statein%vgrs(i,1)*Statein%vgrs(i,1)) & + max(zero, min(Tbd%phy_f2d(i,Model%num_p2d), 30.0)), one) + !*## CCPP ## enddo !*## CCPP ## @@ -1909,6 +1893,30 @@ subroutine GFS_physics_driver & endif !lsm + !! Strictly speaking, this is not required. But when + !! hunting for bit-for-bit differences, updating the + !! subsurface variables in the Sfcprop DDT makes + !! life a lot easier + !if (Model%frac_grid) then + ! do k=1,lsoil + ! do i=1,im + ! if (dry(i)) then + ! Sfcprop%smc(i,k) = smsoil(i,k) + ! Sfcprop%stc(i,k) = stsoil(i,k) + ! Sfcprop%slc(i,k) = slsoil(i,k) + ! endif + ! enddo + ! enddo + !else + ! do k=1,lsoil + ! do i=1,im + ! Sfcprop%smc(i,k) = smsoil(i,k) + ! Sfcprop%stc(i,k) = stsoil(i,k) + ! Sfcprop%slc(i,k) = slsoil(i,k) + ! enddo + ! enddo + !endif + ! if (lprnt) write(0,*)' tseabeficemodel =',Sfcprop%tsfc(ipr),' me=',me & ! &, ' kdt=',kdt,' tsfc32=',tsfc3(ipr,2),' fice=',fice(ipr) & ! &,' stsoil=',stsoil(ipr,:) @@ -1931,8 +1939,9 @@ subroutine GFS_physics_driver & (im, Statein%tgrs(:,1), & Statein%qgrs(:,1,1), cd3(:,2), cdq3(:,2), & Statein%prsl(:,1), wind, & - flag_cice, flag_iter, dqsfc_cice, dtsfc_cice, & - dusfc_cice, dvsfc_cice, & + flag_cice, flag_iter, & + Coupling%dqsfcin_cpl, Coupling%dtsfcin_cpl, & + Coupling%dusfcin_cpl, Coupling%dvsfcin_cpl, & ! --- outputs: qss3(:,2), cmm3(:,2), chh3(:,2), evap3(:,2), hflx3(:,2), & stress3(:,2)) @@ -2103,15 +2112,11 @@ subroutine GFS_physics_driver & ep1d(i) = ep1d3(i,k) Sfcprop%weasd(i) = weasd3(i,k) Sfcprop%snowd(i) = snowd3(i,k) - evap(i) = evap3(i,k) hflx(i) = hflx3(i,k) qss(i) = qss3(i,k) Sfcprop%tsfc(i) = tsfc3(i,k) - Diag%cmm(i) = cmm3(i,k) - Diag%chh(i) = chh3(i,k) - Sfcprop%zorll(i) = zorl3(i,1) Sfcprop%zorlo(i) = zorl3(i,3) @@ -2120,7 +2125,6 @@ subroutine GFS_physics_driver & txo = one - txi evap(i) = txi * evap3(i,2) + txo * evap3(i,3) hflx(i) = txi * hflx3(i,2) + txo * hflx3(i,3) -! Sfcprop%tsfc(i) = txi * tice(i) + txo * tsfc3(i,3) Sfcprop%tsfc(i) = txi * tsfc3(i,2) + txo * tsfc3(i,3) else ! return updated lake ice thickness & concentration to global array if (islmsk(i) == 2) then @@ -2839,10 +2843,10 @@ subroutine GFS_physics_driver & do i=1,im if (Sfcprop%oceanfrac(i) > zero) then ! Ocean only, NO LAKES if (Sfcprop%fice(i) > one - epsln) then ! no open water, thus use results from CICE - Coupling%dusfci_cpl(i) = dusfc_cice(i) - Coupling%dvsfci_cpl(i) = dvsfc_cice(i) - Coupling%dtsfci_cpl(i) = dtsfc_cice(i) - Coupling%dqsfci_cpl(i) = dqsfc_cice(i) + Coupling%dusfci_cpl(i) = Coupling%dusfcin_cpl(i) + Coupling%dvsfci_cpl(i) = Coupling%dvsfcin_cpl(i) + Coupling%dtsfci_cpl(i) = Coupling%dtsfcin_cpl(i) + Coupling%dqsfci_cpl(i) = Coupling%dqsfcin_cpl(i) elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point tem1 = max(Diag%q1(i), 1.e-8) rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1)) @@ -2856,8 +2860,6 @@ subroutine GFS_physics_driver & endif Coupling%dtsfci_cpl(i) = con_cp * rho * hflx3(i,3) ! sensible heat flux over open ocean Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean -! if (lprnt .and. i == ipr) write(0,*)' hflx33=',hflx3(i,3),' evap33=',evap3(i,3), & -! ' con_cp=',con_cp,' rho=',rho,' con_hvap=',con_hvap else ! use results from PBL scheme for 100% open ocean Coupling%dusfci_cpl(i) = dusfc1(i) Coupling%dvsfci_cpl(i) = dvsfc1(i) @@ -5281,7 +5283,7 @@ subroutine GFS_physics_driver & !*## CCPP ## !## CCPP ##* GFS_MP_generic.F90/GFS_MP_generic_post_run Diag%rain(:) = Diag%rainc(:) + frain * rain1(:) ! total rain per timestep - + ! --- get the amount of different precip type for Noah MP ! --- convert from m/dtp to mm/s if (Model%lsm==Model%lsm_noahmp) then @@ -5461,10 +5463,10 @@ subroutine GFS_physics_driver & if (Model%cplflx .or. Model%cplchm) then do i = 1, im - Coupling%rain_cpl(i) = Coupling%rain_cpl(i) & - + Diag%rain(i) * (one-Sfcprop%srflag(i)) - Coupling%snow_cpl(i) = Coupling%snow_cpl(i) & - + Diag%rain(i) * Sfcprop%srflag(i) + Tbd%drain_cpl(i)= Diag%rain(i) * (one-Sfcprop%srflag(i)) + Tbd%dsnow_cpl(i)= Diag%rain(i) * Sfcprop%srflag(i) + Coupling%rain_cpl(i) = Coupling%rain_cpl(i) + Tbd%drain_cpl(i) + Coupling%snow_cpl(i) = Coupling%snow_cpl(i) + Tbd%dsnow_cpl(i) enddo endif @@ -5555,15 +5557,6 @@ subroutine GFS_physics_driver & if (Model%do_sppt) then !--- radiation heating rate Tbd%dtdtr(1:im,:) = Tbd%dtdtr(1:im,:) + dtdtc(1:im,:)*dtf - do i = 1, im - if (t850(i) > 273.16) then -!--- change in change in rain precip - Tbd%drain_cpl(i) = Diag%rain(i) - Tbd%drain_cpl(i) - else -!--- change in change in snow precip - Tbd%dsnow_cpl(i) = Diag%rain(i) - Tbd%dsnow_cpl(i) - endif - enddo endif !*## CCPP ## !## CCPP ##* This block is not in the CCPP since it is not needed in the CCPP. diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 9667ea8e1..d9f420897 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -284,8 +284,8 @@ module GFS_typedefs #endif real (kind=kind_phys), pointer :: q2m (:) => null() !< 2 meter humidity -! -- In/Out for Noah MP - real (kind=kind_phys), pointer :: snowxy (:) => null() ! +! -- In/Out for Noah MP + real (kind=kind_phys), pointer :: snowxy (:) => null() !< real (kind=kind_phys), pointer :: tvxy (:) => null() !< veg temp real (kind=kind_phys), pointer :: tgxy (:) => null() !< ground temp real (kind=kind_phys), pointer :: canicexy(:) => null() !< @@ -312,7 +312,7 @@ module GFS_typedefs real (kind=kind_phys), pointer :: xlaixy (:) => null() !< real (kind=kind_phys), pointer :: taussxy (:) => null() !< real (kind=kind_phys), pointer :: smcwtdxy(:) => null() !< - real (kind=kind_phys), pointer :: deeprechxy(:) => null() !< + real (kind=kind_phys), pointer :: deeprechxy(:)=> null() !< real (kind=kind_phys), pointer :: rechxy (:) => null() !< real (kind=kind_phys), pointer :: snicexy (:,:) => null() !< @@ -1853,11 +1853,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: tsurf_land(:) => null() !< real (kind=kind_phys), pointer :: tsurf_ocean(:) => null() !< real (kind=kind_phys), pointer :: ud_mf(:,:) => null() !< - real (kind=kind_phys), pointer :: ulwsfc_cice(:) => null() !< - real (kind=kind_phys), pointer :: dusfc_cice(:) => null() !< - real (kind=kind_phys), pointer :: dvsfc_cice(:) => null() !< - real (kind=kind_phys), pointer :: dqsfc_cice(:) => null() !< - real (kind=kind_phys), pointer :: dtsfc_cice(:) => null() !< real (kind=kind_phys), pointer :: uustar_ice(:) => null() !< real (kind=kind_phys), pointer :: uustar_land(:) => null() !< real (kind=kind_phys), pointer :: uustar_ocean(:) => null() !< @@ -2582,7 +2577,7 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%ca_turb = clear_val Coupling%ca_shal = clear_val Coupling%ca_rad = clear_val - Coupling%ca_micro = clear_val + Coupling%ca_micro = clear_val Coupling%cape = clear_val Coupling%tconvtend = clear_val Coupling%qconvtend = clear_val @@ -4831,6 +4826,8 @@ subroutine tbd_create (Tbd, IM, Model) if ( Model%isubc_lw == 2 .or. Model%isubc_sw == 2 ) then allocate (Tbd%icsdsw (IM)) allocate (Tbd%icsdlw (IM)) + Tbd%icsdsw = zero + Tbd%icsdlw = zero endif !--- ozone and stratosphere h2o needs @@ -4875,18 +4872,20 @@ subroutine tbd_create (Tbd, IM, Model) Tbd%acvb = clear_val Tbd%acvt = clear_val + if (Model%cplflx .or. Model%cplchm) then + allocate (Tbd%drain_cpl (IM)) + allocate (Tbd%dsnow_cpl (IM)) + Tbd%drain_cpl = clear_val + Tbd%dsnow_cpl = clear_val + endif + if (Model%do_sppt) then allocate (Tbd%dtdtr (IM,Model%levs)) allocate (Tbd%dtotprcp (IM)) allocate (Tbd%dcnvprcp (IM)) - allocate (Tbd%drain_cpl (IM)) - allocate (Tbd%dsnow_cpl (IM)) - Tbd%dtdtr = clear_val Tbd%dtotprcp = clear_val Tbd%dcnvprcp = clear_val - Tbd%drain_cpl = clear_val - Tbd%dsnow_cpl = clear_val endif allocate (Tbd%phy_f2d (IM,Model%ntot2d)) @@ -5425,8 +5424,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%u10max = zero Diag%v10max = zero Diag%spd10max = zero -! Diag%rain = zero -! Diag%rainc = zero + Diag%rain = zero + Diag%rainc = zero Diag%ice = zero Diag%snow = zero Diag%graupel = zero @@ -5909,11 +5908,6 @@ subroutine interstitial_create (Interstitial, IM, Model) allocate (Interstitial%tsurf_land (IM)) allocate (Interstitial%tsurf_ocean (IM)) allocate (Interstitial%ud_mf (IM,Model%levs)) - allocate (Interstitial%ulwsfc_cice (IM)) - allocate (Interstitial%dusfc_cice (IM)) - allocate (Interstitial%dvsfc_cice (IM)) - allocate (Interstitial%dtsfc_cice (IM)) - allocate (Interstitial%dqsfc_cice (IM)) allocate (Interstitial%uustar_ice (IM)) allocate (Interstitial%uustar_land (IM)) allocate (Interstitial%uustar_ocean (IM)) @@ -6436,11 +6430,6 @@ subroutine interstitial_phys_reset (Interstitial, Model) Interstitial%tsurf_land = huge Interstitial%tsurf_ocean = huge Interstitial%ud_mf = clear_val - Interstitial%ulwsfc_cice = clear_val - Interstitial%dusfc_cice = clear_val - Interstitial%dvsfc_cice = clear_val - Interstitial%dtsfc_cice = clear_val - Interstitial%dqsfc_cice = clear_val Interstitial%uustar_ice = huge Interstitial%uustar_land = huge Interstitial%uustar_ocean = huge @@ -6781,11 +6770,6 @@ subroutine interstitial_print(Interstitial, Model, mpirank, omprank, blkno) write (0,*) 'sum(Interstitial%tsurf_land ) = ', sum(Interstitial%tsurf_land ) write (0,*) 'sum(Interstitial%tsurf_ocean ) = ', sum(Interstitial%tsurf_ocean ) write (0,*) 'sum(Interstitial%ud_mf ) = ', sum(Interstitial%ud_mf ) - write (0,*) 'sum(Interstitial%ulwsfc_cice ) = ', sum(Interstitial%ulwsfc_cice ) - write (0,*) 'sum(Interstitial%dusfc_cice ) = ', sum(Interstitial%dusfc_cice ) - write (0,*) 'sum(Interstitial%dvsfc_cice ) = ', sum(Interstitial%dvsfc_cice ) - write (0,*) 'sum(Interstitial%dtsfc_cice ) = ', sum(Interstitial%dtsfc_cice ) - write (0,*) 'sum(Interstitial%dqsfc_cice ) = ', sum(Interstitial%dqsfc_cice ) write (0,*) 'sum(Interstitial%uustar_ice ) = ', sum(Interstitial%uustar_ice ) write (0,*) 'sum(Interstitial%uustar_land ) = ', sum(Interstitial%uustar_land ) write (0,*) 'sum(Interstitial%uustar_ocean ) = ', sum(Interstitial%uustar_ocean ) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.meta b/gfsphysics/GFS_layer/GFS_typedefs.meta index 152a9e8ad..c87a74bf2 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -6579,34 +6579,6 @@ dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys -[dusfc_cice] - standard_name = surface_x_momentum_flux_for_coupling_interstitial - long_name = sfc x momentum flux for coupling interstitial - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys -[dvsfc_cice] - standard_name = surface_y_momentum_flux_for_coupling_interstitial - long_name = sfc y momentum flux for coupling interstitial - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys -[dtsfc_cice] - standard_name = surface_upward_sensible_heat_flux_for_coupling_interstitial - long_name = sfc sensible heat flux for coupling interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys -[dqsfc_cice] - standard_name = surface_upward_latent_heat_flux_for_coupling_interstitial - long_name= surface latent heat flux for coupling interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys [elvmax] standard_name = maximum_subgrid_orography long_name = maximum of subgrid orography @@ -8038,13 +8010,6 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys -[ulwsfc_cice] - standard_name = surface_upwelling_longwave_flux_for_coupling_interstitial - long_name = surface upwelling longwave flux for coupling_interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys [uustar_ocean] standard_name = surface_friction_velocity_over_ocean long_name = surface friction velocity over ocean diff --git a/gfsphysics/physics/GFS_debug.F90 b/gfsphysics/physics/GFS_debug.F90 index 75fa97603..c0b24ca97 100644 --- a/gfsphysics/physics/GFS_debug.F90 +++ b/gfsphysics/physics/GFS_debug.F90 @@ -394,7 +394,12 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Coupling%rain_cpl', Coupling%rain_cpl) call print_var(mpirank,omprank, blkno, 'Coupling%snow_cpl', Coupling%snow_cpl) end if + if (Model%cplwav2atm) then + call print_var(mpirank,omprank, blkno, 'Coupling%zorlwav_cpl' , Coupling%zorlwav_cpl ) + end if if (Model%cplflx) then + call print_var(mpirank,omprank, blkno, 'Coupling%oro_cpl' , Coupling%oro_cpl ) + call print_var(mpirank,omprank, blkno, 'Coupling%slmsk_cpl' , Coupling%slmsk_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%slimskin_cpl', Coupling%slimskin_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%dusfcin_cpl ', Coupling%dusfcin_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%dvsfcin_cpl ', Coupling%dvsfcin_cpl ) @@ -458,11 +463,24 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Coupling%shum_wts', Coupling%shum_wts) end if if (Model%do_skeb) then - call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts) - call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts) + call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts ) + call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts ) end if if (Model%do_sfcperts) then - call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts', Coupling%sfc_wts) + call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts' , Coupling%sfc_wts ) + end if + if (Model%do_ca) then + call print_var(mpirank,omprank, blkno, 'Coupling%tconvtend', Coupling%tconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%qconvtend', Coupling%qconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%uconvtend', Coupling%uconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%vconvtend', Coupling%vconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_out ', Coupling%ca_out ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_deep ', Coupling%ca_deep ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_turb ', Coupling%ca_turb ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_shal ', Coupling%ca_shal ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_rad ', Coupling%ca_rad ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_micro ', Coupling%ca_micro ) + call print_var(mpirank,omprank, blkno, 'Coupling%cape ', Coupling%cape ) end if if(Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then call print_var(mpirank,omprank, blkno, 'Coupling%nwfa2d', Coupling%nwfa2d) diff --git a/gfsphysics/physics/gcycle.F90 b/gfsphysics/physics/gcycle.F90 index e3666c26a..b410aaa9f 100644 --- a/gfsphysics/physics/gcycle.F90 +++ b/gfsphysics/physics/gcycle.F90 @@ -8,7 +8,7 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) GFS_sfcprop_type, GFS_cldprop_type implicit none - integer :: nblks + integer, intent(in) :: nblks type(GFS_control_type), intent(in) :: Model type(GFS_grid_type), intent(in) :: Grid(nblks) type(GFS_sfcprop_type), intent(inout) :: Sfcprop(nblks) @@ -34,7 +34,7 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) TG3FCS (Model%nx*Model%ny), & CNPFCS (Model%nx*Model%ny), & AISFCS (Model%nx*Model%ny), & - F10MFCS(Model%nx*Model%ny), & +! F10MFCS(Model%nx*Model%ny), & VEGFCS (Model%nx*Model%ny), & VETFCS (Model%nx*Model%ny), & SOTFCS (Model%nx*Model%ny), & @@ -103,7 +103,7 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) ZORFCS (len) = Sfcprop(nb)%zorl (ix) TG3FCS (len) = Sfcprop(nb)%tg3 (ix) CNPFCS (len) = Sfcprop(nb)%canopy (ix) - F10MFCS (len) = Sfcprop(nb)%f10m (ix) +! F10MFCS (len) = Sfcprop(nb)%f10m (ix) VEGFCS (len) = Sfcprop(nb)%vfrac (ix) VETFCS (len) = Sfcprop(nb)%vtype (ix) SOTFCS (len) = Sfcprop(nb)%stype (ix) @@ -190,8 +190,8 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) len = len + 1 Sfcprop(nb)%slmsk (ix) = SLIFCS (len) if ( Model%nstf_name(1) > 0 ) then - Sfcprop(nb)%tref(ix) = TSFFCS (len) -! if (Model%nstf_name(2) == 0) then + Sfcprop(nb)%tref(ix) = TSFFCS (len) +! if ( Model%nstf_name(2) == 0 ) then ! dt_warm = (Sfcprop(nb)%xt(ix) + Sfcprop(nb)%xt(ix) ) & ! / Sfcprop(nb)%xz(ix) ! Sfcprop(nb)%tsfco(ix) = Sfcprop(nb)%tref(ix) & @@ -205,7 +205,7 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) Sfcprop(nb)%zorl (ix) = ZORFCS (len) Sfcprop(nb)%tg3 (ix) = TG3FCS (len) Sfcprop(nb)%canopy (ix) = CNPFCS (len) - Sfcprop(nb)%f10m (ix) = F10MFCS (len) +! Sfcprop(nb)%f10m (ix) = F10MFCS (len) Sfcprop(nb)%vfrac (ix) = VEGFCS (len) Sfcprop(nb)%vtype (ix) = VETFCS (len) Sfcprop(nb)%stype (ix) = SOTFCS (len) From 00be176ed5d228063023e11d33f831f438d0b1b1 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Thu, 2 Apr 2020 14:19:44 +0000 Subject: [PATCH 18/21] put gctrt in .no.flxform to avoid debug error for csawmgshoc --- ccpp/physics | 2 +- gfsphysics/physics/cs_conv.F90 | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/ccpp/physics b/ccpp/physics index 84daf0cb3..3d64654af 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 84daf0cb31446e6998dbd1bc40d6bcf6d86e826a +Subproject commit 3d64654aff54e50bdac27f25af9712050b5531d9 diff --git a/gfsphysics/physics/cs_conv.F90 b/gfsphysics/physics/cs_conv.F90 index e824d93ce..8f3b41082 100644 --- a/gfsphysics/physics/cs_conv.F90 +++ b/gfsphysics/physics/cs_conv.F90 @@ -1225,9 +1225,11 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions gcht(i,ctp) = tem * gcht(i,ctp) gcqt(i,ctp) = tem * gcqt(i,ctp) gcit(i,ctp) = tem * gcit(i,ctp) - do n = ntrq,ntr - gctrt(i,n,ctp) = tem * gctrt(i,n,ctp) - enddo + if (.not. flx_form) then + do n = ntrq,ntr + gctrt(i,n,ctp) = tem * gctrt(i,n,ctp) + enddo + end if gcut(i,ctp) = tem * gcut(i,ctp) gcvt(i,ctp) = tem * gcvt(i,ctp) do k=1,kmax From 084b359c2678c0818efabecf91520b6a62c817f3 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Mon, 6 Apr 2020 20:53:53 +0000 Subject: [PATCH 19/21] changes make changing INPUT/cam5_* to cam5_* in iccninterp --- gfsphysics/physics/iccninterp.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gfsphysics/physics/iccninterp.f90 b/gfsphysics/physics/iccninterp.f90 index 66dd344db..d1254692c 100644 --- a/gfsphysics/physics/iccninterp.f90 +++ b/gfsphysics/physics/iccninterp.f90 @@ -33,7 +33,7 @@ SUBROUTINE read_cidata (me, master) end do end do call nf_close(ncid) - call nf_open("INPUT/cam5_4_143_NPCCN_monclimo2.nc", nf_NOWRITE, ncid) + call nf_open("cam5_4_143_NPCCN_monclimo2.nc", nf_NOWRITE, ncid) call nf_inq_varid(ncid, "NPCCN", varid) call nf_get_var(ncid, varid, ccnin) call nf_close(ncid) From 783feb2f92bd76a7bf4194fd03bb478aabffdb75 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Fri, 10 Apr 2020 12:29:35 +0000 Subject: [PATCH 20/21] attached ccpp/physics submodule --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 3d64654af..c1fb9ccb9 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 3d64654aff54e50bdac27f25af9712050b5531d9 +Subproject commit c1fb9ccb98cab43e9006aa3457e1e4a9a393f59d From e82650d66126752b2fa1d01ad33152d361c0bbc5 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Fri, 10 Apr 2020 13:23:32 -0400 Subject: [PATCH 21/21] Revert change to .gitmodules and update submodule pointer for ccpp-physics --- .gitmodules | 2 +- ccpp/physics | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index a4bf59883..d253f6966 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,5 +8,5 @@ branch = master [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/AnningCheng-NOAA/ccpp-physics + url = https://github.com/NCAR/ccpp-physics branch = master diff --git a/ccpp/physics b/ccpp/physics index 3d64654af..b5765fcbf 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 3d64654aff54e50bdac27f25af9712050b5531d9 +Subproject commit b5765fcbf4d6fbb839feb6dc9748470c0f7735df