diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere index 371a29afb..80ce8ce20 160000 --- a/atmos_cubed_sphere +++ b/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit 371a29afbf813357dd93647cac0cbcd44db2ab20 +Subproject commit 80ce8ce200f90985097860c4eb47da1574b87c58 diff --git a/ccpp/physics b/ccpp/physics index efb68b5b9..b5765fcbf 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit efb68b5b948937f256a1a90c2de446b0d9b09e0f +Subproject commit b5765fcbf4d6fbb839feb6dc9748470c0f7735df 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..f1522d7fa --- /dev/null +++ b/ccpp/suites/suite_FV3_GFS_v16_csawmg.xml @@ -0,0 +1,93 @@ + + + + + + + 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 + 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 + GFS_GWD_generic_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 + 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 + + + + diff --git a/gfsphysics/GFS_layer/GFS_diagnostics.F90 b/gfsphysics/GFS_layer/GFS_diagnostics.F90 index 95f7f51e7..7bde5a03a 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_driver.F90 b/gfsphysics/GFS_layer/GFS_driver.F90 index b72225735..c75f7dec9 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -210,10 +210,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 @@ -286,7 +286,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, & @@ -295,7 +295,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, & @@ -1009,7 +1009,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, & @@ -1021,7 +1021,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 98edf87d1..1929f030b 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -4949,32 +4949,30 @@ subroutine GFS_physics_driver & ! write(1000+me,*)' maxwatncb=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt',kdt ! enddo -!## CCPP ##* m_micro.F90/m_micro_run - call m_micro_driver (im, ix, levs, Model%flipv, dtp, Statein%prsl, & - Statein%prsi, Statein%phil, Statein%phii, & - Statein%vvl, clw(1,1,2), QLCN, clw(1,1,1), QICN, & - Radtend%htrlw, Radtend%htrsw, w_upi, cf_upi, & - FRLAND, Diag%HPBL, CNV_MFD, CNV_DQLDT, & -! FRLAND, Diag%HPBL, CNV_MFD, CNV_PRC3, CNV_DQLDT, & - CLCN, Stateout%gu0, Stateout%gv0, Diag%dusfc, & - Diag%dvsfc, dusfc1, dvsfc1, dusfc1, dvsfc1, & - CNV_FICE, CNV_NDROP, CNV_NICE, Stateout%gq0(1,1,1), & - Stateout%gq0(1,1,ntcw), & - Stateout%gq0(1,1,ntiw), Stateout%gt0, rain1, & - Diag%sr, Stateout%gq0(1,1,ntlnc), & - Stateout%gq0(1,1,ntinc), Model%fprcp, qrn, & - qsnw, qgl, ncpr, ncps, ncgl, & - Tbd%phy_f3d(1,1,1), kbot, & - 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, & - skip_macro, lprnt, & -! skip_macro, cn_prc, cn_snr, lprnt, & -! ipr, kdt, Grid%xlat, Grid%xlon) - Model%mg_alf, Model%mg_qcmin, Model%pdfflag, & - ipr, kdt, Grid%xlat, Grid%xlon, rhc) -!*## CCPP ## + call m_micro_driver (im, ix, levs, Model%flipv, dtp, Statein%prsl, & + Statein%prsi, Statein%phil, Statein%phii, & + Statein%vvl, clw(1,1,2), QLCN, clw(1,1,1), QICN, & + Radtend%htrlw, Radtend%htrsw, w_upi, cf_upi, & + FRLAND, Diag%HPBL, CNV_MFD, CNV_DQLDT, & +! FRLAND, Diag%HPBL, CNV_MFD, CNV_PRC3, CNV_DQLDT, & + CLCN, Stateout%gu0, Stateout%gv0, Diag%dusfc, & + Diag%dvsfc, dusfc1, dvsfc1, dusfc1, dvsfc1, & + CNV_FICE, CNV_NDROP, CNV_NICE, Stateout%gq0(1,1,1), & + Stateout%gq0(1,1,ntcw), & + Stateout%gq0(1,1,ntiw), Stateout%gt0, rain1, & + Diag%sr, Stateout%gq0(1,1,ntlnc), & + Stateout%gq0(1,1,ntinc), Model%fprcp, qrn, & + qsnw, qgl, ncpr, ncps, ncgl, & + Tbd%phy_f3d(1,1,1), kbot, & + 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, & + Tbd%in_nm, Tbd%ccn_nm, Model%iccn, & + skip_macro, lprnt, & +! skip_macro, cn_prc, cn_snr, lprnt, & +! ipr, kdt, Grid%xlat, Grid%xlon) + Model%mg_alf, Model%mg_qcmin, Model%pdfflag, & + ipr, kdt, Grid%xlat, Grid%xlon, rhc) ! do k=1,levs ! write(1000+me,*)' maxwatnca=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt=',kdt ! enddo diff --git a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 index d14eeaac3..dad425a73 100644 --- a/gfsphysics/GFS_layer/GFS_radiation_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_radiation_driver.F90 @@ -1566,7 +1566,8 @@ subroutine GFS_radiation_driver & !check print *,' in grrad : calling setaer ' !## CCPP ##* GFS_rrtmg_pre.F90/GFS_rrtmg_pre_run 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 @@ -2058,12 +2059,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/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 00bae3897..d9f420897 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -612,7 +612,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 @@ -1066,7 +1066,7 @@ module GFS_typedefs real(kind=kind_phys) :: julian !< julian day using midnight of January 1 of forecast year as initial epoch integer :: yearlen !< length of the current forecast year in days ! - 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 @@ -2735,8 +2735,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) @@ -3104,7 +3104,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- coupling parameters cplflx, cplwav, cplwav2atm, 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 @@ -3320,17 +3320,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) @@ -3340,6 +3330,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 @@ -3365,7 +3364,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 @@ -4229,7 +4228,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, & @@ -4434,7 +4433,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 @@ -4784,7 +4782,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)) @@ -4794,7 +4792,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 bc2344250..c87a74bf2 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.meta +++ b/gfsphysics/GFS_layer/GFS_typedefs.meta @@ -2113,9 +2113,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 @@ -3858,9 +3858,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/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..6a0eadb99 100644 --- a/gfsphysics/physics/aerinterp.f90 +++ b/gfsphysics/physics/aerinterp.f90 @@ -1,165 +1,187 @@ 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) + 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) + + 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 +219,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 +237,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 +260,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 +282,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 +302,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 +326,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/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 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) diff --git a/gfsphysics/physics/m_micro_driver.F90 b/gfsphysics/physics/m_micro_driver.F90 index 26a04d96a..8801a05c2 100644 --- a/gfsphysics/physics/m_micro_driver.F90 +++ b/gfsphysics/physics/m_micro_driver.F90 @@ -13,7 +13,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 & &, lprnt, alf_fac, qc_min, pdfflag & &, ipr, kdt, xlat, xlon, rhc_i) @@ -60,7 +60,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) :: & @@ -300,7 +301,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 @@ -361,7 +362,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 @@ -531,7 +532,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 @@ -586,7 +587,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 @@ -769,12 +770,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 @@ -837,19 +833,21 @@ 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 (.not. iccn) then + if (iccn.ne.1) then CDNC_NUC(I,k) = npccninr8(k) INC_NUC (I,k) = naair8(k) endif @@ -984,7 +982,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) @@ -1133,11 +1131,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 @@ -1351,7 +1345,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) @@ -1487,7 +1481,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 281802878..ec44523e6 100644 --- a/gfsphysics/physics/micro_mg2_0.F90 +++ b/gfsphysics/physics/micro_mg2_0.F90 @@ -393,7 +393,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, & @@ -464,7 +464,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 @@ -1113,7 +1114,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) @@ -1152,7 +1153,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 @@ -1167,7 +1168,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 f27aa1896..a9df06c6c 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,8 +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 ! Inputs for aerosol activation @@ -1441,7 +1441,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 +1495,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,11 +1510,13 @@ 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 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 @@ -2844,7 +2846,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)- & diff --git a/gfsphysics/physics/radiation_aerosols.f b/gfsphysics/physics/radiation_aerosols.f index 37364b8be..45a909ca8 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 : ntrcaerm + ! 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 /= ntrcaerm ) 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) ------------------------! ! ======================================================================= !..........................................! diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index e9358ec91..d4ebf3fec 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -1161,6 +1161,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & endif enddo enddo +! print *,'in gfs_post, get tisfc=',maxval(ti), minval(ti) endif ! sea ice skin temperature