diff --git a/.github/workflows/GCC.yml b/.github/workflows/GCC.yml index 44d827db08..7e63c3c609 100644 --- a/.github/workflows/GCC.yml +++ b/.github/workflows/GCC.yml @@ -120,7 +120,7 @@ jobs: export CC=mpicc export CXX=mpicxx export FC=mpif90 - cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON ${{ matrix.cmake_opts }} ${{ env.gcov_cmake }} + cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON -DFV3=ON ${{ matrix.cmake_opts }} ${{ env.gcov_cmake }} make -j2 - name: run-tests diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 1137d4b881..8d4a185f5e 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -119,7 +119,7 @@ jobs: export CC=mpicc export CXX=mpicxx export FC=mpif90 - cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON ${{ matrix.cmake_opts }} -DENABLE_FV3ATM_DOCS=ON ${{ env.gcov_cmake }} + cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON -DFV3=ON ${{ matrix.cmake_opts }} -DENABLE_FV3ATM_DOCS=ON ${{ env.gcov_cmake }} make doxygen_doc - name: upload-docs diff --git a/.gitmodules b/.gitmodules index 19bde8db17..c8ea487ea5 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,8 +8,8 @@ branch = main [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/grantfirl/ccpp-physics - branch = feature/rad-fix_gjf + url = https://github.com/dustinswales/ccpp-physics + branch = feature/rad-fix_djs [submodule "upp"] path = upp url = https://github.com/NOAA-EMC/UPP diff --git a/ccpp/config/ccpp_prebuild_config.py b/ccpp/config/ccpp_prebuild_config.py index 8fa2460d54..968c48a772 100755 --- a/ccpp/config/ccpp_prebuild_config.py +++ b/ccpp/config/ccpp_prebuild_config.py @@ -20,6 +20,7 @@ 'physics/physics/Radiation/RRTMG/radsw_param.f', 'physics/physics/Radiation/RRTMG/radlw_param.f', 'physics/physics/photochem/module_ozphys.F90', + 'physics/physics/MP/TEMPO/TEMPO/module_mp_tempo_params.F90', 'physics/physics/photochem/module_h2ophys.F90', 'physics/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.F90', '../ccpp/data/CCPP_typedefs.F90', @@ -46,6 +47,10 @@ 'module_ozphys' : '', 'ty_ozphys' : '', }, + 'module_mp_tempo_params' : { + 'module_mp_tempo_params' : '', + 'ty_tempo_cfg' : '', + }, 'module_h2ophys' : { 'module_h2ophys' : '', 'ty_h2ophys' : '', @@ -162,8 +167,9 @@ 'physics/physics/photochem/module_h2ophys.F90', 'physics/physics/photochem/module_ozphys.F90', 'physics/physics/MP/Ferrier_Aligo/mp_fer_hires.F90', - 'physics/physics/MP/GFDL/gfdl_cloud_microphys.F90', + 'physics/physics/MP/GFDL/v1_2019/gfdl_cloud_microphys.F90', 'physics/physics/MP/GFDL/fv_sat_adj.F90', + 'physics/physics/MP/GFDL/v3_2022/gfdl_cloud_microphys_v3.F90', 'physics/physics/MP/Morrison_Gettelman/m_micro.F90', 'physics/physics/MP/Morrison_Gettelman/m_micro_pre.F90', 'physics/physics/MP/Morrison_Gettelman/m_micro_post.F90', @@ -171,6 +177,9 @@ 'physics/physics/MP/Thompson/mp_thompson_pre.F90', 'physics/physics/MP/Thompson/mp_thompson.F90', 'physics/physics/MP/Thompson/mp_thompson_post.F90', + 'physics/physics/MP/TEMPO/mp_tempo_pre.F90', + 'physics/physics/MP/TEMPO/mp_tempo.F90', + 'physics/physics/MP/TEMPO/mp_tempo_post.F90', 'physics/physics/MP/Zhao_Carr/zhaocarr_gscond.f', 'physics/physics/MP/Zhao_Carr/zhaocarr_precpd.f', 'physics/physics/PBL/HEDMF/hedmf.f', diff --git a/ccpp/data/CCPP_typedefs.F90 b/ccpp/data/CCPP_typedefs.F90 index 8b24554ccf..d9aaebf3b2 100644 --- a/ccpp/data/CCPP_typedefs.F90 +++ b/ccpp/data/CCPP_typedefs.F90 @@ -820,7 +820,7 @@ subroutine gfs_interstitial_create (Interstitial, IM, Model) ! ! Allocate arrays that are conditional on physics choices if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson & - .or. Model%imp_physics == Model%imp_physics_nssl & + .or. Model%imp_physics == Model%imp_physics_tempo .or. Model%imp_physics == Model%imp_physics_nssl & ) then allocate (Interstitial%graupelmp (IM)) allocate (Interstitial%icemp (IM)) @@ -908,7 +908,8 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model) ! perform aerosol convective transport and PBL diffusion Interstitial%trans_aero = Model%cplchm .and. Model%trans_trac - if (Model%imp_physics == Model%imp_physics_thompson) then + if (Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) then if (Model%ltaerosol) then Interstitial%nvdiff = 12 else if (Model%mraerosol) then @@ -961,7 +962,8 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model) if (Model%imp_physics == Model%imp_physics_wsm6) then Interstitial%ntcwx = 2 Interstitial%ntiwx = 3 - elseif (Model%imp_physics == Model%imp_physics_thompson) then + elseif (Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) then Interstitial%ntcwx = 2 Interstitial%ntiwx = 3 Interstitial%ntrwx = 4 @@ -999,7 +1001,8 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model) endif elseif (Model%imp_physics == Model%imp_physics_gfdl) then Interstitial%nvdiff = 7 - elseif (Model%imp_physics == Model%imp_physics_thompson) then + elseif (Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) then if (Model%ltaerosol) then Interstitial%nvdiff = 12 else if (Model%mraerosol) then @@ -1421,7 +1424,7 @@ subroutine gfs_interstitial_phys_reset (Interstitial, Model) ! ! Reset fields that are conditional on physics choices if (Model%imp_physics == Model%imp_physics_gfdl .or. Model%imp_physics == Model%imp_physics_thompson & - .or. Model%imp_physics == Model%imp_physics_nssl & + .or. Model%imp_physics == Model%imp_physics_tempo .or. Model%imp_physics == Model%imp_physics_nssl & ) then Interstitial%graupelmp = clear_val Interstitial%icemp = clear_val diff --git a/ccpp/data/CCPP_typedefs.meta b/ccpp/data/CCPP_typedefs.meta index c4de321ca9..6e6f8e90c1 100644 --- a/ccpp/data/CCPP_typedefs.meta +++ b/ccpp/data/CCPP_typedefs.meta @@ -1128,7 +1128,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme) + active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) [gwdcu] standard_name = tendency_of_x_wind_due_to_convective_gravity_wave_drag long_name = zonal wind tendency due to convective gravity wave drag @@ -1220,7 +1220,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme) + active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) [dry] standard_name = flag_nonzero_land_surface_fraction long_name = flag indicating presence of some land surface area fraction @@ -1711,7 +1711,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme) + active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) [rainp] standard_name = tendency_of_rain_water_mixing_ratio_due_to_microphysics long_name = tendency of rain water mixing ratio due to microphysics @@ -1967,7 +1967,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme) + active = (control_for_microphysics_scheme == identifier_for_gfdl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_nssl_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) [snowmt] standard_name = surface_snow_melt long_name = snow melt during timestep diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90 index 8d90054332..a5febdbb28 100644 --- a/ccpp/data/GFS_typedefs.F90 +++ b/ccpp/data/GFS_typedefs.F90 @@ -11,10 +11,12 @@ module GFS_typedefs con_omega, con_rerth, con_psat, karman, rainmin,& con_c, con_plnk, con_boltz, con_solr_2008, & con_solr_2002, con_thgni, con_1ovg, con_rgas, & - con_avgd, con_amd, con_amw + con_avgd, con_amd, con_amw, con_one, con_p001, & + con_secinday use module_radsw_parameters, only: topfsw_type, sfcfsw_type use module_radlw_parameters, only: topflw_type, sfcflw_type + use module_mp_tempo_params, only: ty_tempo_cfg use module_ozphys, only: ty_ozphys use module_h2ophys, only: ty_h2ophys use land_iau_mod, only: land_iau_external_data_type, land_iau_control_type, & @@ -38,7 +40,8 @@ module GFS_typedefs integer, parameter :: dfi_radar_max_intervals = 4 !< Number of radar-derived temperature tendency and/or convection suppression intervals. Do not change. real(kind=kind_phys), parameter :: limit_unspecified = 1e12 !< special constant for "namelist value was not provided" in radar-derived temperature tendency limit range - + + integer, parameter :: physics_no_tracer = -99 !> \section arg_table_GFS_typedefs !! \htmlinclude GFS_typedefs.html @@ -169,6 +172,11 @@ module GFS_typedefs real (kind=kind_phys), pointer :: vvl (:,:) => null() !< layer mean vertical velocity in pa/sec real (kind=kind_phys), pointer :: tgrs (:,:) => null() !< model layer mean temperature in k real (kind=kind_phys), pointer :: qgrs (:,:,:) => null() !< layer mean tracer concentration +!3D-SA-TKE + real (kind=kind_phys), pointer :: def_1 (:,:) => null() !< deformation + real (kind=kind_phys), pointer :: def_2 (:,:) => null() !< deformation + real (kind=kind_phys), pointer :: def_3 (:,:) => null() !< deformation +!3D-SA-TKE-end ! dissipation estimate real (kind=kind_phys), pointer :: diss_est(:,:) => null() !< model layer mean temperature in k ! soil state variables - for soil SPPT - sfc-perts, mgehne @@ -944,6 +952,7 @@ module GFS_typedefs integer :: imp_physics !< choice of microphysics scheme integer :: imp_physics_gfdl = 11 !< choice of GFDL microphysics scheme integer :: imp_physics_thompson = 8 !< choice of Thompson microphysics scheme + integer :: imp_physics_tempo = 88 !< choice of TEMPO microphysics scheme integer :: imp_physics_wsm6 = 6 !< choice of WSMG microphysics scheme integer :: imp_physics_zhao_carr = 99 !< choice of Zhao-Carr microphysics scheme integer :: imp_physics_zhao_carr_pdf = 98 !< choice of Zhao-Carr microphysics scheme with PDF clouds @@ -1028,6 +1037,7 @@ module GFS_typedefs !--- Thompson's microphysical parameters logical :: ltaerosol !< flag for aerosol version logical :: mraerosol !< flag for merra2_aerosol_aware + logical :: lthailaware !< flag for TEMPO hail-aware logical :: lradar !< flag for radar reflectivity real(kind=kind_phys) :: nsfullradar_diag!< seconds between resetting radar reflectivity calculation real(kind=kind_phys) :: ttendlim !< temperature tendency limiter per time step in K/s @@ -1036,6 +1046,7 @@ module GFS_typedefs real(kind=kind_phys) :: dt_inner !< time step for the inner loop in s logical :: sedi_semi !< flag for semi Lagrangian sedi of rain integer :: decfl !< deformed CFL factor + type(ty_tempo_cfg) :: tempo_cfg !< Thompson MP configuration information. logical :: thompson_mp_is_init=.false. !< Local scheme initialization flag !--- GFDL microphysical paramters @@ -1156,6 +1167,7 @@ module GFS_typedefs logical :: do_ugwp_v1 !< flag for version 1 ugwp GWD logical :: do_ugwp_v1_orog_only !< flag for version 1 ugwp GWD (orographic drag only) logical :: do_ugwp_v1_w_gsldrag !< flag for version 1 ugwp with OGWD of GSL + logical :: do_ngw_ec !< flag for ecmwf ngw logical :: mstrat !< flag for moorthi approach for stratus logical :: moist_adj !< flag for moist convective adjustment logical :: cscnv !< flag for Chikira-Sugiyama convection @@ -1181,6 +1193,7 @@ module GFS_typedefs logical :: shinhong !< flag for scale-aware Shinhong vertical turbulent mixing scheme logical :: do_ysu !< flag for YSU turbulent mixing scheme logical :: dspheat !< flag for tke dissipative heating + logical :: sa3dtke !< flag for scale-aware 3D tke scheme logical :: hurr_pbl !< flag for hurricane-specific options in PBL scheme logical :: lheatstrg !< flag for canopy heat storage parameterization logical :: lseaspray !< flag for sea spray parameterization @@ -1377,6 +1390,7 @@ module GFS_typedefs real(kind=kind_phys) :: elmx !< maximum allowed dissipation mixing length in boundary layer mass flux scheme integer :: sfc_rlm !< choice of near surface mixing length in boundary layer mass flux scheme integer :: tc_pbl !< control for TC applications in the PBL scheme + integer :: use_lpt !< control for using Liquid Potential Temp for TC applications in the GFSPBL scheme !--- parameters for canopy heat storage (CHS) parameterization real(kind=kind_phys) :: h0facu !< CHS factor for sensible heat flux in unstable surface layer @@ -1653,6 +1667,8 @@ module GFS_typedefs real(kind=kind_phys) :: rhcmax ! maximum critical relative humidity, replaces rhc_max in physcons.F90 real(kind=kind_phys) :: huge !< huge fill value +!--- AQM Canopy + logical :: do_canopy !< control flag for aqm canopy effects !--- lightning threat and diagsnostics logical :: lightning_threat !< report lightning threat indices @@ -1798,6 +1814,11 @@ module GFS_typedefs real (kind=kind_phys), pointer :: hpbl (:) => null() !< Planetary boundary layer height real (kind=kind_phys), pointer :: ud_mf (:,:) => null() !< updraft mass flux +!-- Diagnostic variable that passes to dyn_core (SA-3D-TKE) + real (kind=kind_phys), pointer :: dku3d_h (:,:) => null() !< Horizontal eddy diffusitivity for momentum + real (kind=kind_phys), pointer :: dku3d_e (:,:) => null() !< Eddy diffusitivity for momentum for tke + + !--- dynamical forcing variables for Grell-Freitas convection real (kind=kind_phys), pointer :: forcet (:,:) => null() !< real (kind=kind_phys), pointer :: forceq (:,:) => null() !< @@ -2193,6 +2214,18 @@ module GFS_typedefs ! Diagnostics for coupled air quality model real (kind=kind_phys), pointer :: aod (:) => null() !< instantaneous aerosol optical depth ( n/a ) +!IVAI + ! Diagnostics for coupled air quality model + real (kind=kind_phys), pointer :: coszens(:) => null() ! Cosine SZA for photolysis + real (kind=kind_phys), pointer :: jo3o1d(:) => null() ! instantaneous O3O1D photolysis rate + real (kind=kind_phys), pointer :: jno2 (:) => null() ! instantaneous NO2 photolysis rate + real (kind=kind_phys), pointer :: claie(:) => null() ! Leaf Area Index ECCC + real (kind=kind_phys), pointer :: cfch (:) => null() ! Forest Canopy Height + real (kind=kind_phys), pointer :: cfrt (:) => null() ! Forest Fraction + real (kind=kind_phys), pointer :: cclu (:) => null() ! Clumping Index + real (kind=kind_phys), pointer :: cpopu(:) => null() ! Population density +!IVAI + ! Auxiliary output arrays for debugging real (kind=kind_phys), pointer :: aux2d(:,:) => null() !< auxiliary 2d arrays in output (for debugging) real (kind=kind_phys), pointer :: aux3d(:,:,:)=> null() !< auxiliary 2d arrays in output (for debugging) @@ -2208,6 +2241,10 @@ module GFS_typedefs real (kind=kind_phys), pointer :: do3_dt_temp(:,:) => null() real (kind=kind_phys), pointer :: do3_dt_ohoz(:,:) => null() + !--- NRL WV photochemistry diagnostics + real (kind=kind_phys), pointer :: dqv_dt_prd(:, :) => null() + real (kind=kind_phys), pointer :: dqv_dt_qvmx(:, :) => null() + contains procedure :: create => diag_create procedure :: rad_zero => diag_rad_zero @@ -2273,11 +2310,21 @@ subroutine statein_create (Statein, Model) allocate (Statein%wgrs (IM,Model%levs)) endif allocate (Statein%qgrs (IM,Model%levs,Model%ntrac)) +!3D-SA-TKE + allocate (Statein%def_1 (IM,Model%levs)) + allocate (Statein%def_2 (IM,Model%levs)) + allocate (Statein%def_3 (IM,Model%levs)) +!3D-SA-TKE-end Statein%qgrs = clear_val Statein%pgr = clear_val Statein%ugrs = clear_val Statein%vgrs = clear_val +!3D-SA-TKE + Statein%def_1 = clear_val + Statein%def_2 = clear_val + Statein%def_3 = clear_val +!3D-SA-TKE-end if(Model%lightning_threat) then Statein%wgrs = clear_val @@ -3268,7 +3315,9 @@ subroutine coupling_create (Coupling, Model) endif !--- needed for Thompson's aerosol option - if(Model%imp_physics == Model%imp_physics_thompson .and. (Model%ltaerosol .or. Model%mraerosol)) then + if((Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) .and. & + (Model%ltaerosol .or. Model%mraerosol)) then allocate (Coupling%nwfa2d (IM)) allocate (Coupling%nifa2d (IM)) Coupling%nwfa2d = clear_val @@ -3328,7 +3377,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- modules use physcons, only: con_rerth, con_pi use mersenne_twister, only: random_setseed, random_number - use parse_tracers, only: get_tracer_index ! implicit none @@ -3438,7 +3486,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iems = 0 !< 1.0 => Noah lsm !< 2.0 => Noah MP and RUC lsms integer :: iaer = 1 !< default aerosol effect in sw only - integer :: iaermdl = 0 !< default tropospheric aerosol model scheme flag + integer :: iaermdl = 1 !< default tropospheric aerosol model scheme flag !< 0: seasonal global distributed OPAC aerosol climatology !< 1: monthly global distributed GOCART aerosol climatology !< 2: GOCART prognostic aerosol model @@ -3595,6 +3643,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- Thompson microphysical parameters logical :: ltaerosol = .false. !< flag for aerosol version logical :: mraerosol = .false. !< flag for merra2_aerosol_aware + logical :: lthailaware = .false. !< flag for TEMPO hail-aware logical :: lradar = .false. !< flag for radar reflectivity real(kind=kind_phys) :: nsfullradar_diag = -999.0 !< seconds between resetting radar reflectivity calculation, set to <0 for every time step real(kind=kind_phys) :: ttendlim = -999.0 !< temperature tendency limiter, set to <0 to deactivate @@ -3708,6 +3757,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: do_ugwp_v1_orog_only = .false. !< flag for version 1 ugwp GWD (orographic drag only) logical :: do_ugwp_v1_w_gsldrag = .false. !< flag for version 1 ugwp GWD (orographic drag only) !--- vay-2018 + logical :: do_ngw_ec = .false. !< flag for ecmwf ngwd algorithm logical :: ldiag_ugwp = .false. !< flag for UGWP diag fields logical :: ugwp_seq_update = .false. !< flag for updating winds between UGWP steps logical :: do_ugwp = .false. !< flag do UGWP+RF @@ -3736,6 +3786,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: shinhong = .false. !< flag for scale-aware Shinhong vertical turbulent mixing scheme logical :: do_ysu = .false. !< flag for YSU vertical turbulent mixing scheme logical :: dspheat = .false. !< flag for tke dissipative heating + logical :: sa3dtke = .false. !< flag for scale-aware 3D tke scheme logical :: hurr_pbl = .false. !< flag for hurricane-specific options in PBL scheme logical :: lheatstrg = .false. !< flag for canopy heat storage parameterization logical :: lseaspray = .false. !< flag for sea spray parameterization @@ -3918,6 +3969,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & real(kind=kind_phys) :: elmx = 300. !< maximum allowed dissipation mixing length in boundary layer mass flux scheme integer :: sfc_rlm = 0 !< choice of near surface mixing length in boundary layer mass flux scheme integer :: tc_pbl = 0 !< control for TC applications in the PBL scheme + integer :: use_lpt = 0 !< control for using Liquid Potential Temp for TC applications in the GFSPBL scheme !--- parameters for canopy heat storage (CHS) parameterization real(kind=kind_phys) :: h0facu = 0.25 @@ -4062,6 +4114,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! logical :: land_iau_do_stcsmc_adjustment = .false. ! real(kind=kind_phys) :: land_iau_min_T_increment = 0.0001 +!--- switch for aqm canopy effects + logical :: do_canopy = .false. !< flag for canopy option + !--- END NAMELIST VARIABLES NAMELIST /gfs_physics_nml/ & @@ -4104,8 +4159,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & mg_do_graupel, mg_do_hail, mg_nccons, mg_nicons, mg_ngcons, & mg_ncnst, mg_ninst, mg_ngnst, sed_supersat, do_sb_physics, & mg_alf, mg_qcmin, mg_do_ice_gmao, mg_do_liq_liu, & - ltaerosol, lradar, nsfullradar_diag, lrefres, ttendlim, & - ext_diag_thompson, dt_inner, lgfdlmprad, & + ltaerosol, lthailaware, lradar, nsfullradar_diag, lrefres, & + ttendlim, ext_diag_thompson, dt_inner, lgfdlmprad, & sedi_semi, decfl, & nssl_cccn, nssl_alphah, nssl_alphahl, & nssl_alphar, nssl_ehw0, nssl_ehlw0, & @@ -4144,7 +4199,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & gwd_opt, do_ugwp_v0, do_ugwp_v0_orog_only, & do_ugwp_v0_nst_only, & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & - do_gwd_opt_psl, & + do_gwd_opt_psl, do_ngw_ec, & do_ugwp_v1, do_ugwp_v1_orog_only, do_ugwp_v1_w_gsldrag, & ugwp_seq_update, var_ric, coef_ric_l, coef_ric_s, hurr_pbl, & do_myjsfc, do_myjpbl, & @@ -4162,6 +4217,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & do_spp, n_var_spp, & lndp_type, n_var_lndp, lndp_each_step, & pert_mp,pert_clds,pert_radtend, & + !--- Scale-aware 3D TKE scheme + sa3dtke, & !--- Rayleigh friction prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, & ! --- Ferrier-Aligo @@ -4184,7 +4241,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & diag_flux, diag_log, & ! vertical diffusion xkzm_m, xkzm_h, xkzm_s, xkzminv, moninq_fac, dspfac, & - bl_upfr, bl_dnfr, rlmx, elmx, sfc_rlm, tc_pbl, & + bl_upfr, bl_dnfr, rlmx, elmx, sfc_rlm, tc_pbl, use_lpt, & !--- canopy heat storage parameterization h0facu, h0facs, & !--- cellular automata @@ -4219,6 +4276,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- (DFI) time ranges with radar-prescribed microphysics tendencies ! and (maybe) convection suppression fh_dfi_radar, radar_tten_limits, do_cap_suppress, & + ! aqm canopy option + do_canopy, & !--- GSL lightning threat indices lightning_threat !, & ! !--- land_iau_nml @@ -4576,8 +4635,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ialb = ialb Model%iems = iems Model%iaer = iaer + Model%iaermdl = iaer/1000 Model%iaerclm = iaerclm - if (iaer/1000 == 1 .or. Model%iccn == 2) then + if (iaer/1000 == 1 .or. Model%iccn == 2 .or. Model%iaermdl ==6 ) then Model%iaerclm = .true. ntrcaer = ntrcaerm else if (iaer/1000 == 2) then @@ -4587,7 +4647,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & endif Model%lalw1bd = lalw1bd Model%iaerflg = iaerflg - Model%iaermdl = iaermdl Model%aeros_file = aeros_file Model%solar_file = solar_file Model%semis_file = semis_file @@ -4765,6 +4824,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- Thompson MP parameters Model%ltaerosol = ltaerosol Model%mraerosol = mraerosol + Model%lthailaware = lthailaware if (Model%ltaerosol .and. Model%mraerosol) then write(0,*) 'Logic error: Only one Thompson aerosol option can be true, either ltaerosol or mraerosol)' stop @@ -4780,6 +4840,16 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & endif Model%sedi_semi = sedi_semi Model%decfl = decfl + +!--- TEMPO MP parameters + ! DJS to Anders: Maybe we put more of these nml options into the TEMPO configuration type? + Model%tempo_cfg%aerosol_aware = (ltaerosol .or. mraerosol) + Model%tempo_cfg%hail_aware = lthailaware + if (Model%ltaerosol .and. Model%mraerosol) then + write(0,*) 'Logic error: Only one TEMPO aerosol option can be true, either ltaerosol or mraerosol)' + stop + end if + !--- F-A MP parameters Model%rhgrd = rhgrd Model%spec_adv = spec_adv @@ -4883,7 +4953,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%exticeden = exticeden if (Model%exticeden .and. & (Model%imp_physics /= Model%imp_physics_gfdl .and. Model%imp_physics /= Model%imp_physics_thompson .and. & - Model%imp_physics /= Model%imp_physics_nssl )) then + Model%imp_physics /= Model%imp_physics_nssl .and. Model%imp_physics /= Model%imp_physics_tempo)) then !see GFS_MP_generic_post.F90; exticeden is only compatible with GFDL, !Thompson, or NSSL MP print *,' Using exticeden = T is only valid when using GFDL, Thompson, or NSSL microphysics.' @@ -5072,6 +5142,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%do_gsl_drag_ss = do_gsl_drag_ss Model%do_gsl_drag_tofd = do_gsl_drag_tofd Model%do_gwd_opt_psl = do_gwd_opt_psl + Model%do_ngw_ec = do_ngw_ec Model%do_ugwp_v1 = do_ugwp_v1 Model%do_ugwp_v1_orog_only = do_ugwp_v1_orog_only Model%do_ugwp_v1_w_gsldrag = do_ugwp_v1_w_gsldrag @@ -5141,6 +5212,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%diag_flux = diag_flux !--- flux method in 2-m diagnostics (for stable conditions) Model%diag_log = diag_log +!--- SA-3D-TKE option + Model%sa3dtke = sa3dtke !--- vertical diffusion Model%xkzm_m = xkzm_m @@ -5155,6 +5228,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%elmx = elmx Model%sfc_rlm = sfc_rlm Model%tc_pbl = tc_pbl + Model%use_lpt = use_lpt !--- canopy heat storage parametrization Model%h0facu = h0facu @@ -5178,6 +5252,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%do_spp = do_spp Model%n_var_spp = n_var_spp +!--- aqm canopy effects in physics + Model%do_canopy = do_canopy + if (Model%lndp_type/=0) then allocate (Model%lndp_var_list(Model%n_var_lndp)) allocate (Model%lndp_prt_list(Model%n_var_lndp)) @@ -5245,51 +5322,50 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%tracer_names(:) = tracer_names(:) Model%ntqv = 1 #ifdef MULTI_GASES - Model%nto = get_tracer_index(Model%tracer_names, 'spo', Model%me, Model%master, Model%debug) - Model%nto2 = get_tracer_index(Model%tracer_names, 'spo2', Model%me, Model%master, Model%debug) - Model%ntoz = get_tracer_index(Model%tracer_names, 'spo3', Model%me, Model%master, Model%debug) + Model%nto = get_physics_tracer_index('spo', Model) + Model%nto2 = get_physics_tracer_index('spo2', Model) + Model%ntoz = get_physics_tracer_index('spo3', Model) #else - Model%ntoz = get_tracer_index(Model%tracer_names, 'o3mr', Model%me, Model%master, Model%debug) + Model%ntoz = get_physics_tracer_index('o3mr', Model) if( Model%ntoz <= 0 ) & - Model%ntoz = get_tracer_index(Model%tracer_names, 'spo3', Model%me, Model%master, Model%debug) + Model%ntoz = get_physics_tracer_index('spo3', Model) #endif - Model%ntcw = get_tracer_index(Model%tracer_names, 'liq_wat', Model%me, Model%master, Model%debug) - Model%ntiw = get_tracer_index(Model%tracer_names, 'ice_wat', Model%me, Model%master, Model%debug) - Model%ntrw = get_tracer_index(Model%tracer_names, 'rainwat', Model%me, Model%master, Model%debug) - Model%ntsw = get_tracer_index(Model%tracer_names, 'snowwat', Model%me, Model%master, Model%debug) - Model%ntgl = get_tracer_index(Model%tracer_names, 'graupel', Model%me, Model%master, Model%debug) - Model%nthl = get_tracer_index(Model%tracer_names, 'hailwat', Model%me, Model%master, Model%debug) - Model%ntclamt = get_tracer_index(Model%tracer_names, 'cld_amt', Model%me, Model%master, Model%debug) - Model%ntlnc = get_tracer_index(Model%tracer_names, 'water_nc', Model%me, Model%master, Model%debug) - Model%ntinc = get_tracer_index(Model%tracer_names, 'ice_nc', Model%me, Model%master, Model%debug) - Model%ntrnc = get_tracer_index(Model%tracer_names, 'rain_nc', Model%me, Model%master, Model%debug) - Model%ntsnc = get_tracer_index(Model%tracer_names, 'snow_nc', Model%me, Model%master, Model%debug) - Model%ntgnc = get_tracer_index(Model%tracer_names, 'graupel_nc', Model%me, Model%master, Model%debug) - Model%nthnc = get_tracer_index(Model%tracer_names, 'hail_nc', Model%me, Model%master, Model%debug) - Model%ntccn = get_tracer_index(Model%tracer_names, 'ccn_nc', Model%me, Model%master, Model%debug) - Model%ntccna = get_tracer_index(Model%tracer_names, 'ccna_nc', Model%me, Model%master, Model%debug) - Model%ntgv = get_tracer_index(Model%tracer_names, 'graupel_vol',Model%me, Model%master, Model%debug) - Model%nthv = get_tracer_index(Model%tracer_names, 'hail_vol', Model%me, Model%master, Model%debug) - Model%ntrz = get_tracer_index(Model%tracer_names, 'rain_ref', Model%me, Model%master, Model%debug) - Model%ntgz = get_tracer_index(Model%tracer_names, 'graupel_ref',Model%me, Model%master, Model%debug) - Model%nthz = get_tracer_index(Model%tracer_names, 'hail_ref', Model%me, Model%master, Model%debug) - Model%ntke = get_tracer_index(Model%tracer_names, 'sgs_tke', Model%me, Model%master, Model%debug) - Model%ntsigma = get_tracer_index(Model%tracer_names, 'sigmab', Model%me, Model%master, Model%debug) - Model%ntomega = get_tracer_index(Model%tracer_names, 'omegab', Model%me, Model%master, Model%debug) - Model%nqrimef = get_tracer_index(Model%tracer_names, 'q_rimef', Model%me, Model%master, Model%debug) - Model%ntwa = get_tracer_index(Model%tracer_names, 'liq_aero', Model%me, Model%master, Model%debug) - Model%ntia = get_tracer_index(Model%tracer_names, 'ice_aero', Model%me, Model%master, Model%debug) + Model%ntcw = get_physics_tracer_index('liq_wat', Model) + Model%ntiw = get_physics_tracer_index('ice_wat', Model) + Model%ntrw = get_physics_tracer_index('rainwat', Model) + Model%ntsw = get_physics_tracer_index('snowwat', Model) + Model%ntgl = get_physics_tracer_index('graupel', Model) + Model%nthl = get_physics_tracer_index('hailwat', Model) + Model%ntclamt = get_physics_tracer_index('cld_amt', Model) + Model%ntlnc = get_physics_tracer_index('water_nc', Model) + Model%ntinc = get_physics_tracer_index('ice_nc', Model) + Model%ntrnc = get_physics_tracer_index('rain_nc', Model) + Model%ntsnc = get_physics_tracer_index('snow_nc', Model) + Model%ntgnc = get_physics_tracer_index('graupel_nc', Model) + Model%nthnc = get_physics_tracer_index('hail_nc', Model) + Model%ntccn = get_physics_tracer_index('ccn_nc', Model) + Model%ntccna = get_physics_tracer_index('ccna_nc', Model) + Model%ntgv = get_physics_tracer_index('graupel_vol', Model) + Model%nthv = get_physics_tracer_index('hail_vol', Model) + Model%ntrz = get_physics_tracer_index('rain_ref', Model) + Model%ntgz = get_physics_tracer_index('graupel_ref', Model) + Model%nthz = get_physics_tracer_index('hail_ref', Model) + Model%ntke = get_physics_tracer_index('sgs_tke', Model) + Model%ntsigma = get_physics_tracer_index('sigmab', Model) + Model%ntomega = get_physics_tracer_index('omegab', Model) + Model%nqrimef = get_physics_tracer_index('q_rimef', Model) + Model%ntwa = get_physics_tracer_index('liq_aero', Model) + Model%ntia = get_physics_tracer_index('ice_aero', Model) if (Model%cpl_fire) then - Model%ntfsmoke = get_tracer_index(Model%tracer_names, 'fsmoke', Model%me, Model%master, Model%debug) + Model%ntfsmoke = get_physics_tracer_index('fsmoke', Model) endif if (Model%rrfs_sd) then - Model%ntsmoke = get_tracer_index(Model%tracer_names, 'smoke', Model%me, Model%master, Model%debug) - Model%ntdust = get_tracer_index(Model%tracer_names, 'dust', Model%me, Model%master, Model%debug) - Model%ntcoarsepm = get_tracer_index(Model%tracer_names, 'coarsepm', Model%me, Model%master, Model%debug) + Model%ntsmoke = get_physics_tracer_index('smoke', Model) + Model%ntdust = get_physics_tracer_index('dust', Model) + Model%ntcoarsepm = get_physics_tracer_index('coarsepm', Model) endif !--- initialize parameters for atmospheric chemistry tracers call Model%init_chemistry(tracer_types) - !--- setup aerosol scavenging factors call Model%init_scavenging(fscav_aero) @@ -5340,26 +5416,27 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Last index of outermost dimension of dtend Model%ndtend = 0 allocate (Model%dtidx(Model%ntracp100,Model%nprocess)) - Model%dtidx = -99 + Model%dtidx = physics_no_tracer if(Model%ntchm>0) then - Model%ntdu1 = get_tracer_index(Model%tracer_names, 'dust1', Model%me, Model%master, Model%debug) - Model%ntdu2 = get_tracer_index(Model%tracer_names, 'dust2', Model%me, Model%master, Model%debug) - Model%ntdu3 = get_tracer_index(Model%tracer_names, 'dust3', Model%me, Model%master, Model%debug) - Model%ntdu4 = get_tracer_index(Model%tracer_names, 'dust4', Model%me, Model%master, Model%debug) - Model%ntdu5 = get_tracer_index(Model%tracer_names, 'dust5', Model%me, Model%master, Model%debug) - Model%ntss1 = get_tracer_index(Model%tracer_names, 'seas1', Model%me, Model%master, Model%debug) - Model%ntss2 = get_tracer_index(Model%tracer_names, 'seas2', Model%me, Model%master, Model%debug) - Model%ntss3 = get_tracer_index(Model%tracer_names, 'seas3', Model%me, Model%master, Model%debug) - Model%ntss4 = get_tracer_index(Model%tracer_names, 'seas4', Model%me, Model%master, Model%debug) - Model%ntss5 = get_tracer_index(Model%tracer_names, 'seas5', Model%me, Model%master, Model%debug) - Model%ntsu = get_tracer_index(Model%tracer_names, 'so4', Model%me, Model%master, Model%debug) - Model%ntbcb = get_tracer_index(Model%tracer_names, 'bc1', Model%me, Model%master, Model%debug) - Model%ntbcl = get_tracer_index(Model%tracer_names, 'bc2', Model%me, Model%master, Model%debug) - Model%ntocb = get_tracer_index(Model%tracer_names, 'oc1', Model%me, Model%master, Model%debug) - Model%ntocl = get_tracer_index(Model%tracer_names, 'oc2', Model%me, Model%master, Model%debug) + Model%ntdu1 = get_physics_tracer_index('dust1', Model) + Model%ntdu2 = get_physics_tracer_index('dust2', Model) + Model%ntdu3 = get_physics_tracer_index('dust3', Model) + Model%ntdu4 = get_physics_tracer_index('dust4', Model) + Model%ntdu5 = get_physics_tracer_index('dust5', Model) + Model%ntss1 = get_physics_tracer_index('seas1', Model) + Model%ntss2 = get_physics_tracer_index('seas2', Model) + Model%ntss3 = get_physics_tracer_index('seas3', Model) + Model%ntss4 = get_physics_tracer_index('seas4', Model) + Model%ntss5 = get_physics_tracer_index('seas5', Model) + Model%ntsu = get_physics_tracer_index('so4', Model) + Model%ntbcb = get_physics_tracer_index('bc1', Model) + Model%ntbcl = get_physics_tracer_index('bc2', Model) + Model%ntocb = get_physics_tracer_index('oc1', Model) + Model%ntocl = get_physics_tracer_index('oc2', Model) end if + ! Lake & fractional grid safety checks if(Model%me==Model%master) then if(Model%lkm>0 .and. Model%frac_grid) then @@ -5421,17 +5498,18 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Need better descriptions of these. call label_dtend_tracer(Model,100+Model%ntchm+Model%ntchs-1,'pp10','pp10 concentration','kg kg-1 s-1') - itrac=get_tracer_index(Model%tracer_names, 'DMS', Model%me, Model%master, Model%debug) + itrac=get_physics_tracer_index('DMS', Model) if(itrac>0) then call label_dtend_tracer(Model,100+itrac,'DMS','DMS concentration','kg kg-1 s-1') endif - itrac=get_tracer_index(Model%tracer_names, 'msa', Model%me, Model%master, Model%debug) + itrac=get_physics_tracer_index('msa', Model) if(itrac>0) then call label_dtend_tracer(Model,100+itrac,'msa','msa concentration','kg kg-1 s-1') endif endif endif + call label_dtend_tracer(Model,Model%index_of_temperature,'temp','temperature','K s-1') call label_dtend_tracer(Model,Model%index_of_x_wind,'u','x wind','m s-2') call label_dtend_tracer(Model,Model%index_of_y_wind,'v','y wind','m s-2') @@ -5464,7 +5542,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & call label_dtend_tracer(Model,100+Model%ntia,'ice_aero','number concentration of ice-friendly aerosols','kg-1 s-1') call label_dtend_tracer(Model,100+Model%nto,'o_ion','oxygen ion concentration','kg kg-1 s-1') call label_dtend_tracer(Model,100+Model%nto2,'o2','oxygen concentration','kg kg-1 s-1') - call label_dtend_cause(Model,Model%index_of_process_pbl,'pbl','tendency due to PBL') call label_dtend_cause(Model,Model%index_of_process_dcnv,'deepcnv','tendency due to deep convection') call label_dtend_cause(Model,Model%index_of_process_scnv,'shalcnv','tendency due to shallow convection') @@ -5514,7 +5591,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & call fill_dtidx(Model,dtend_select,Model%index_of_y_wind,Model%index_of_process_physics) call fill_dtidx(Model,dtend_select,Model%index_of_x_wind,Model%index_of_process_non_physics) call fill_dtidx(Model,dtend_select,Model%index_of_y_wind,Model%index_of_process_non_physics) - if(qdiag3d) then call fill_dtidx(Model,dtend_select,100+Model%ntqv,Model%index_of_process_scnv,have_scnv) call fill_dtidx(Model,dtend_select,100+Model%ntqv,Model%index_of_process_dcnv,have_dcnv) @@ -5557,6 +5633,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & call fill_dtidx(Model,dtend_select,100+Model%ntoz,Model%index_of_process_physics,.true.) call fill_dtidx(Model,dtend_select,100+Model%ntoz,Model%index_of_process_non_physics,.true.) + call fill_dtidx(Model, dtend_select, 100+Model%ntqv,Model%index_of_process_prod_loss, h2o_phys) + call fill_dtidx(Model, dtend_select, 100+Model%ntqv,Model%index_of_process_ozmix, h2o_phys) + call fill_dtidx(Model, dtend_select, 100+Model%ntqv,Model%index_of_process_photochem, h2o_phys) + if(.not.Model%do_mynnedmf .and. .not. Model%satmedmf) then call fill_dtidx(Model,dtend_select,100+Model%ntqv,Model%index_of_process_pbl,have_pbl) call fill_dtidx(Model,dtend_select,100+Model%ntcw,Model%index_of_process_pbl,have_pbl) @@ -5578,7 +5658,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & enddo endif end if - IF ( Model%imp_physics == Model%imp_physics_nssl2mccn ) THEN ! recognize this option for compatibility Model%imp_physics = Model%imp_physics_nssl nssl_ccn_on = .true. @@ -5598,8 +5677,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ENDIF stop ENDIF - Model%ntccn = -99 - Model%ntccna = -99 + Model%ntccn = physics_no_tracer + Model%ntccna = physics_no_tracer ELSEIF ( Model%ntccn < 1 ) THEN if (Model%me == Model%master) then write(*,*) 'NSSL micro: error! CCN is ON but ntccn < 1. Must have ccn_nc in field_table if nssl_ccn_on=T' @@ -5721,7 +5800,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%levh2o = 1 Model%h2o_coeff = 1 end if - !--- quantities to be used to derive phy_f*d totals Model%nshoc_2d = nshoc_2d Model%nshoc_3d = nshoc_3d @@ -5811,8 +5889,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- BEGIN CODE FROM COMPNS_PHYSICS !--- shoc scheme if (do_shoc) then - if (Model%imp_physics == Model%imp_physics_thompson) then - print *,'SHOC is not currently compatible with Thompson MP -- shutting down' + if ((Model%imp_physics == Model%imp_physics_thompson) .or. & + (Model%imp_physics == Model%imp_physics_tempo)) then + print *,'SHOC is not currently compatible with Thompson/TEMPO MP -- shutting down' stop endif Model%nshoc_3d = 3 @@ -6172,7 +6251,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ' num_p2d =',Model%num_p2d - elseif (Model%imp_physics == Model%imp_physics_thompson) then !Thompson microphysics + elseif (Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) then !Thompson/TEMPO microphysics Model%npdf3d = 0 Model%num_p3d = 3 Model%num_p2d = 1 @@ -6189,9 +6269,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & print *,' Thompson MP requires effr_in to be set to .true. - job aborted' stop end if - if (Model%me == Model%master) print *,' Using Thompson double moment microphysics', & + if (Model%me == Model%master) print *,' Using Thompson/TEMPO double moment microphysics', & ' ltaerosol = ',Model%ltaerosol, & ' mraerosol = ',Model%mraerosol, & + ' lthailaware = ',Model%lthailaware, & ' ttendlim =',Model%ttendlim, & ' ext_diag_thompson =',Model%ext_diag_thompson, & ' dt_inner =',Model%dt_inner, & @@ -6364,7 +6445,6 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Model%upd_slc = land_iau_upd_slc ! Model%do_stcsmc_adjustment = land_iau_do_stcsmc_adjustment ! Model%min_T_increment = land_iau_min_T_increment - call Model%print () end subroutine control_initialize @@ -6441,8 +6521,6 @@ subroutine control_chemistry_initialize(Model, tracer_types) !--- prognostic and diagnostic chemistry tracers. !--- Each tracer set is assumed to be contiguous. - use parse_tracers, only: NO_TRACER - !--- interface variables class(GFS_control_type) :: Model integer, intent(in) :: tracer_types(:) @@ -6454,11 +6532,11 @@ subroutine control_chemistry_initialize(Model, tracer_types) Model%nchem = 0 Model%ndvel = 0 Model%ntchm = 0 - Model%ntchs = NO_TRACER - Model%ntche = NO_TRACER + Model%ntchs = physics_no_tracer + Model%ntche = physics_no_tracer Model%ndchm = 0 - Model%ndchs = NO_TRACER - Model%ndche = NO_TRACER + Model%ndchs = physics_no_tracer + Model%ndche = physics_no_tracer if (Model%rrfs_sd) then Model%nchem = 3 @@ -6491,8 +6569,6 @@ end subroutine control_chemistry_initialize !---------------------------- subroutine control_scavenging_initialize(Model, fscav) - use parse_tracers, only: get_tracer_index - !--- interface variables class(GFS_control_type) :: Model character(len=*), intent(in) :: fscav(:) @@ -6526,13 +6602,11 @@ subroutine control_scavenging_initialize(Model, fscav) if (j > 1) then read(fscav(i)(j+1:), *, iostat=ios) tem if (ios /= 0) cycle - n = get_tracer_index(Model%tracer_names, adjustl(fscav(i)(:j-1)), Model%me, Model%master, Model%debug) & - - Model%ntchs + 1 + n = get_physics_tracer_index(adjustl(fscav(i)(:j-1)), Model) - Model%ntchs + 1 if (n > 0) Model%fscav(n) = tem endif enddo endif - end subroutine control_scavenging_initialize @@ -6746,10 +6820,12 @@ subroutine control_print(Model) print *, ' wminco : ', Model%wminco print *, ' ' endif - if (Model%imp_physics == Model%imp_physics_wsm6 .or. Model%imp_physics == Model%imp_physics_thompson) then + if ((Model%imp_physics == Model%imp_physics_wsm6) .or. (Model%imp_physics == Model%imp_physics_thompson) .or. & + (Model%imp_physics == Model%imp_physics_tempo)) then print *, ' Thompson microphysical parameters' print *, ' ltaerosol : ', Model%ltaerosol print *, ' mraerosol : ', Model%mraerosol + print *, ' lthailaware : ', Model%lthailaware print *, ' lradar : ', Model%lradar print *, ' nsfullradar_diag : ', Model%nsfullradar_diag print *, ' lrefres : ', Model%lrefres @@ -6904,6 +6980,7 @@ subroutine control_print(Model) print *, ' shinhong : ', Model%shinhong print *, ' do_ysu : ', Model%do_ysu print *, ' dspheat : ', Model%dspheat + print *, ' sa3dtke : ', Model%sa3dtke print *, ' lheatstrg : ', Model%lheatstrg print *, ' lseaspray : ', Model%lseaspray print *, ' cnvcld : ', Model%cnvcld @@ -6944,6 +7021,7 @@ subroutine control_print(Model) print *, ' do_gsl_drag_tofd : ', Model%do_gsl_drag_tofd print *, ' do_gwd_opt_psl : ', Model%do_gwd_opt_psl print *, ' do_ugwp_v1 : ', Model%do_ugwp_v1 + print *, ' do_ngw_ec : ', Model%do_ngw_ec print *, ' do_ugwp_v1_orog_only : ', Model%do_ugwp_v1_orog_only print *, ' do_ugwp_v1_w_gsldrag : ', Model%do_ugwp_v1_w_gsldrag print *, ' hurr_pbl : ', Model%hurr_pbl @@ -7000,6 +7078,7 @@ subroutine control_print(Model) print *, ' elmx : ', Model%elmx print *, ' sfc_rlm : ', Model%sfc_rlm print *, ' tc_pbl : ', Model%tc_pbl + print *, ' use_lpt : ', Model%use_lpt print *, ' ' print *, 'parameters for canopy heat storage parametrization' print *, ' h0facu : ', Model%h0facu @@ -7134,6 +7213,7 @@ subroutine control_print(Model) print *, ' first_time_step : ', Model%first_time_step print *, ' restart : ', Model%restart print *, ' lsm_cold_start : ', Model%lsm_cold_start + print *, ' do_canopy : ', Model%do_canopy print *, ' ' print *, 'lightning threat indexes' print *, ' lightning_threat : ', Model%lightning_threat @@ -7352,6 +7432,12 @@ subroutine tbd_create (Tbd, Model) allocate (Tbd%hpbl (IM)) Tbd%hpbl = clear_val +! Allocate horizontal component of dku for dyn_core (SA-3D-TKE) + allocate (Tbd%dku3d_h (IM,Model%levs)) + Tbd%dku3d_h = clear_val + allocate (Tbd%dku3d_e (IM,Model%levs)) + Tbd%dku3d_e = clear_val + if (Model%imfdeepcnv == Model%imfdeepcnv_gf .or. Model%imfdeepcnv == Model%imfdeepcnv_ntiedtke .or. Model%imfdeepcnv == Model%imfdeepcnv_samf .or. Model%imfshalcnv == Model%imfshalcnv_samf .or. Model%imfdeepcnv == Model%imfdeepcnv_c3 .or. Model%imfshalcnv == Model%imfshalcnv_c3) then allocate (Tbd%prevsq(IM, Model%levs)) Tbd%prevsq = clear_val @@ -7747,7 +7833,6 @@ end subroutine label_dtend_cause ! GFS_diag%create !---------------- subroutine diag_create (Diag, Model) - use parse_tracers, only: get_tracer_index class(GFS_diag_type) :: Diag type(GFS_control_type), intent(in) :: Model integer :: IM @@ -7908,6 +7993,10 @@ subroutine diag_create (Diag, Model) allocate (Diag%do3_dt_temp(IM, Model%levs)) allocate (Diag%do3_dt_ohoz(IM, Model%levs)) endif + if (Model%h2o_phys) then + allocate (Diag%dqv_dt_prd( IM, Model%levs)) + allocate (Diag%dqv_dt_qvmx(IM, Model%levs)) + end if endif ! UGWP @@ -8053,6 +8142,42 @@ subroutine diag_create (Diag, Model) Diag%aod = zero end if +!IVAI: + ! Air quality diagnostics + ! -- initialize diagnostic variables + if (Model%cplaqm) then + +!IVAI: photdiag arrays + allocate (Diag%coszens(IM)) + Diag%coszens= zero + + allocate (Diag%jo3o1d(IM)) + Diag%jo3o1d = zero + + allocate (Diag%jno2(IM)) + Diag%jno2 = zero + +!IVAI: canopy arrays read via aqm_emis_read + if (Model%do_canopy) then + allocate (Diag%claie(IM)) + Diag%claie = zero + + allocate (Diag%cfch (IM)) + Diag%cfch = zero + + allocate (Diag%cfrt (IM)) + Diag%cfrt = zero + + allocate (Diag%cclu (IM)) + Diag%cclu = zero + + allocate (Diag%cpopu (IM)) + Diag%cpopu = zero + end if! (Model%do_canopy) + + end if ! (Model%cplaqm) +!IVAI + ! Auxiliary arrays in output for debugging if (Model%naux2d>0) then allocate (Diag%aux2d(IM,Model%naux2d)) @@ -8246,6 +8371,10 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%do3_dt_temp = zero Diag%do3_dt_ohoz = zero endif + if (Model%h2o_phys) then + Diag%dqv_dt_prd = zero + Diag%dqv_dt_qvmx = zero + end if endif ! @@ -8348,5 +8477,22 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) endif end subroutine diag_phys_zero + + function get_physics_tracer_index (name, Model) + !This function uses the FMS version of get_tracer_index, but changes the missing tracer index to the value used throughout the physics code, rather than the one used in FMS + use tracer_manager_mod, only: get_tracer_index, NO_TRACER + use field_manager_mod, only: MODEL_ATMOS + + character(len=*), intent(in) :: name + type(GFS_control_type), intent(in) :: Model + + !--- local variables + integer :: get_physics_tracer_index + + get_physics_tracer_index = get_tracer_index(MODEL_ATMOS, name, verbose = (Model%me == Model%master) .and. Model%debug) + + if (get_physics_tracer_index == NO_TRACER) get_physics_tracer_index = physics_no_tracer + + end function get_physics_tracer_index end module GFS_typedefs diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta index 04562b7db7..096d7c4cc0 100644 --- a/ccpp/data/GFS_typedefs.meta +++ b/ccpp/data/GFS_typedefs.meta @@ -119,6 +119,48 @@ type = real kind = kind_phys active = (do_lightning_threat_index_calculations) +[def_1] + standard_name = square_of_vertical_shear_due_to_dynamics + long_name = square of vertical shear calculated from dynamics + units = m2 s-2 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys +[def_1(:,1)] + standard_name = square_of_vertical_shear_due_to_dynamics_at_surface_adjacent_layer + long_name = square of vertical shear calculated from dynamics at lowest model layer + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys +[def_2] + standard_name = square_of_horizontal_shear_due_to_dynamics + long_name = square of horizontal shear calculated from dynamics + units = m2 s-2 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys +[def_2(:,1)] + standard_name = square_of_horizontal_shear_due_to_dynamics_at_surface_adjacent_layer + long_name = square of horizontal shear calculated from dynamics at lowest model layer + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys +[def_3] + standard_name = horizontal_transfer_rate_of_tke_due_to_dynamics + long_name = rate of horizontal TKE transfer and pressure correlation calculated from dynamics + units = m2 s-3 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys +[def_3(:,1)] + standard_name = horizontal_transfer_rate_tke_due_to_dynamics_at_surface_adjacent_layer + long_name = rate of horizontal TKE transfer and pressure correlation calculated from dynamics at lowest model layer + units = m2 s-3 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys [vvl] standard_name = lagrangian_tendency_of_air_pressure long_name = layer mean vertical velocity @@ -352,7 +394,7 @@ [qgrs(:,:,index_of_updraft_velocity_in_tracer_concentration_array)] standard_name = prognostic_updraft_velocity_in_convection long_name = convective updraft velocity - units = frac + units = Pa s-1 dimensions = (horizontal_dimension,vertical_layer_dimension) type = real kind = kind_phys @@ -642,7 +684,7 @@ [gq0(:,:,index_of_updraft_velocity_in_tracer_concentration_array)] standard_name = updraft_velocity_updated_by_physics long_name = convective updraft area fraction updated by physics - units = frac + units = Pa s-1 dimensions = (horizontal_dimension,vertical_layer_dimension) type = real kind = kind_phys @@ -3177,7 +3219,7 @@ dimensions = (horizontal_dimension) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .and. (flag_for_aerosol_physics .or. do_merra2_aerosol_awareness)) + active = ((control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) .and. (flag_for_aerosol_physics .or. do_merra2_aerosol_awareness)) [nifa2d] standard_name = tendency_of_nonhygroscopic_ice_nucleating_aerosols_at_surface_adjacent_layer long_name = instantaneous ice-friendly sfc aerosol source @@ -3185,7 +3227,7 @@ dimensions = (horizontal_dimension) type = real kind = kind_phys - active = (control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .and. (flag_for_aerosol_physics .or. do_merra2_aerosol_awareness)) + active = ((control_for_microphysics_scheme == identifier_for_thompson_microphysics_scheme .or. control_for_microphysics_scheme == identifier_for_tempo_microphysics_scheme) .and. (flag_for_aerosol_physics .or. do_merra2_aerosol_awareness)) [ebu_smoke] standard_name = ebu_smoke long_name = buffer of vertical fire emission @@ -4359,6 +4401,12 @@ units = flag dimensions = () type = integer +[imp_physics_tempo] + standard_name = identifier_for_tempo_microphysics_scheme + long_name = choice of TEMPO microphysics scheme + units = flag + dimensions = () + type = integer [imp_physics_wsm6] standard_name = identifier_for_wsm6_microphysics_scheme long_name = choice of WSM6 microphysics scheme @@ -4883,6 +4931,12 @@ units = flag dimensions = () type = logical +[lthailaware] + standard_name = flag_for_hail_physics + long_name = flag for hail physics + units = flag + dimensions = () + type = logical [mraerosol] standard_name = do_merra2_aerosol_awareness long_name = flag for merra2 aerosol-aware physics for example the thompson microphysics @@ -4933,6 +4987,12 @@ units = count dimensions = () type = integer +[tempo_cfg] + standard_name = configuration_for_TEMPO_microphysics + long_name = configuration information for TEMPO microphysics + units = mixed + dimensions = () + type = ty_tempo_cfg [thompson_mp_is_init] standard_name = flag_for_thompson_mp_scheme_initialization long_name = flag carrying scheme initialization status @@ -5525,6 +5585,12 @@ units = flag dimensions = () type = logical +[sa3dtke] + standard_name = do_scale_aware_3d_tke + long_name = flag for scale-aware 3d tke scheme + units = flag + dimensions = () + type = logical [hurr_pbl] standard_name = flag_hurricane_PBL long_name = flag for hurricane-specific options in PBL scheme @@ -6112,6 +6178,12 @@ units = none dimensions = () type = integer +[use_lpt] + standard_name = control_for_using_LPT_for_TC_applications_in_the_PBL_scheme + long_name = control for using LPT in TC applications in the PBL scheme + units = none + dimensions = () + type = integer [h0facu] standard_name = multiplicative_tuning_parameter_for_reduced_surface_heat_fluxes_due_to_canopy_heat_storage long_name = canopy heat storage factor for sensible heat flux in unstable surface layer @@ -7663,6 +7735,12 @@ dimensions = () type = real kind = kind_phys +[do_ngw_ec] + standard_name = flag_for_ngw_ec + long_name = flag to activate ecmwf ngwd + units = flag + dimensions = () + type = logical [do_ugwp_v1] standard_name = flag_for_ugwp_version_1 long_name = flag to activate ver 1 CIRES UGWP @@ -7687,6 +7765,12 @@ units = flag dimensions = () type = logical +[do_canopy] + standard_name = flag_for_canopy_option + long_name = flag for in-canopy eddy diffusivity adjustment option + units = flag + dimensions = () + type = logical [clm_lake_depth_default] standard_name = default_lake_depth_in_clm_lake_model long_name = default lake depth in clm lake model @@ -7718,7 +7802,6 @@ units = count dimensions = () type = integer - ######################################################################## [ccpp-table-properties] name = GFS_grid_type @@ -8005,6 +8088,20 @@ dimensions = (horizontal_dimension) type = real kind = kind_phys +[dku3d_h] + standard_name = horizontal_atmosphere_momentum_diffusivity_for_dynamics + long_name = horizontal atmospheric momentum diffusivity for dynamics + units = m2 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys +[dku3d_e] + standard_name = horizontal_atmosphere_tke_diffusivity_for_dynamics + long_name = horizontal atmospheric tke diffusivity for dynamics + units = m2 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys [ud_mf] standard_name = instantaneous_atmosphere_updraft_convective_mass_flux long_name = (updraft mass flux) * delt @@ -8643,6 +8740,46 @@ dimensions = (horizontal_dimension,number_of_diagnostics_variables_for_radiation) type = real kind = kind_phys +[claie] + standard_name = canopy_leaf_area_index + long_name = canopy leaf area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + active = (flag_for_air_quality_coupling .and. flag_for_canopy_option) +[cfch] + standard_name = canopy_forest_height + long_name = canopy forest height + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + active = (flag_for_air_quality_coupling .and. flag_for_canopy_option) +[cfrt] + standard_name = canopy_forest_fraction + long_name = canopy forest fraction + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + active = (flag_for_air_quality_coupling .and. flag_for_canopy_option) +[cclu] + standard_name = canopy_clumping_index + long_name = canopy clumping index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + active = (flag_for_air_quality_coupling .and. flag_for_canopy_option) +[cpopu] + standard_name = canopy_population_density + long_name = population density used for canopy correction + units = km-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + active = (flag_for_air_quality_coupling .and. flag_for_canopy_option) [topfsw] standard_name = sw_fluxes_top_atmosphere long_name = sw radiation fluxes at toa @@ -9473,6 +9610,22 @@ type = real kind = kind_phys active = (flag_for_diagnostics_3D .and. flag_for_nrl_2015_ozone_scheme) +[dqv_dt_prd] + standard_name = water_vapor_tendency_due_to_production_and_loss_rate + long_name = water_vapor tendency due to production and loss rate + units = kg kg-1 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + active = (flag_for_diagnostics_3D .and. flag_for_stratospheric_water_vapor_physics) +[dqv_dt_qvmx] + standard_name = water_vapor_tendency_due_to_water_vapor_mixing_ratio + long_name = water_vapor tendency due to water_vapor mixing ratio + units = kg kg-1 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + active = (flag_for_diagnostics_3D .and. flag_for_stratospheric_water_vapor_physics) [refl_10cm] standard_name = radar_reflectivity_10cm long_name = instantaneous refl_10cm @@ -10180,8 +10333,9 @@ relative_path = ../physics/physics/ dependencies = hooks/machine.F,hooks/physcons.F90 dependencies = Radiation/RRTMG/radlw_param.f,Radiation/RRTMG/radsw_param.f + dependencies = MP/TEMPO/TEMPO/module_mp_tempo_params.F90 dependencies = photochem/module_ozphys.F90,photochem/module_h2ophys.F90 - dependencies = SFC_Models/Land/Noahmp/lnd_iau_mod.F90,MP/GFDL/GFDL_parse_tracers.F90 + dependencies = SFC_Models/Land/Noahmp/lnd_iau_mod.F90 [ccpp-arg-table] name = GFS_typedefs @@ -10560,3 +10714,24 @@ dimensions = () type = real kind = kind_phys +[con_one] + standard_name = constant_one + long_name = mathematical constant of one + units = 1 + dimensions = () + type = real + kind = kind_phys +[con_p001] + standard_name = constant_one_hundredth + long_name = mathematical constant for one hundredth + units = 1 + dimensions = () + type = real + kind = kind_phys +[con_secinday] + standard_name = seconds_in_a_day + long_name = number of seconds in a day + units = s + dimensions = () + type = real + kind = kind_phys diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90 index 9b2757703a..e09d9f7760 100644 --- a/ccpp/driver/GFS_diagnostics.F90 +++ b/ccpp/driver/GFS_diagnostics.F90 @@ -117,7 +117,6 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ! ExtDiag%data%var3(:,:) [real*8 ] pointer to 3D data [=> null() for a 2D field] ! !---------------------------------------------------------------------------------------------! - use parse_tracers, only: get_tracer_index implicit none ! ! --- interface variables @@ -778,6 +777,97 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop endif endif +!IVAI +!--- air quality diagnostics --- + if (Model%cplaqm) then + +! IVAI: photdiag fields + if (associated(IntDiag%coszens)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'COSZENS' + ExtDiag(idx)%desc = 'Cosine Solar Zenith Angle for Photolysis' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%coszens(:) + endif + + if (associated(IntDiag%jo3o1d)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'JO3O1D' + ExtDiag(idx)%desc = 'photolysis rate O3 for canopy correction' + ExtDiag(idx)%unit = 'min-1' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%jo3o1d(:) + endif + + if (associated(IntDiag%jno2)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'JNO2' + ExtDiag(idx)%desc = 'photolysis rate NO2 for canopy correction' + ExtDiag(idx)%unit = 'min-1' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%jno2(:) + endif + +!IVAI: canopy arrays read via aqm_emis_read + if (Model%do_canopy) then + if (associated(IntDiag%claie)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'CLAIE' + ExtDiag(idx)%desc = 'Leaf Area Index ECCC' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%claie(:) + endif + + if (associated(IntDiag%cfch)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'CFCH' + ExtDiag(idx)%desc = 'Forest Canopy Height' + ExtDiag(idx)%unit = 'm' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%cfch(:) + endif + + if (associated(IntDiag%cfrt)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'CFRT' + ExtDiag(idx)%desc = 'Forest Canopy Fraction' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%cfrt(:) + endif + + if (associated(IntDiag%cclu)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'CCLU' + ExtDiag(idx)%desc = 'Canopy Clumping Index' + ExtDiag(idx)%unit = 'numerical' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%cclu(:) + endif + + if (associated(IntDiag%cpopu)) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'CPOPU' + ExtDiag(idx)%desc = 'Population Density for canopy correction' + ExtDiag(idx)%unit = 'km-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + ExtDiag(idx)%data%var2 => IntDiag%cpopu(:) + endif + endif ! (Model%do_canopy) + + end if ! (Model%cplaqm) +!IVAI + ! ! !--- accumulated diagnostics --- @@ -3741,7 +3831,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop endif endif extended_smoke_dust_diagnostics - + idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'ebu_smoke' @@ -3852,7 +3942,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ! Cloud effective radii from Microphysics if (Model%imp_physics == Model%imp_physics_thompson .or. Model%imp_physics == Model%imp_physics_fer_hires .or. & - Model%imp_physics == Model%imp_physics_nssl ) then + Model%imp_physics == Model%imp_physics_nssl .or. Model%imp_physics == Model%imp_physics_tempo ) then idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'cleffr' diff --git a/ccpp/driver/GFS_restart.F90 b/ccpp/driver/GFS_restart.F90 index 2af5373956..fcca0afd97 100644 --- a/ccpp/driver/GFS_restart.F90 +++ b/ccpp/driver/GFS_restart.F90 @@ -140,7 +140,8 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, num2d = num2d + 2 endif ! Thompson aerosol-aware - if (Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then + if ((Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) .and. (Model%ltaerosol)) then num2d = num2d + 2 endif if (Model%do_cap_suppress .and. Model%num_dfi_radar>0) then @@ -421,7 +422,8 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, Restart(idx)%data%var2 => Sfcprop%rainncprv(:) endif ! Thompson aerosol-aware - if (Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then + if ((Model%imp_physics == Model%imp_physics_thompson .or. & + Model%imp_physics == Model%imp_physics_tempo) .and. Model%ltaerosol) then idx = idx + 1 Restart(idx)%name = 'thompson_2d_nwfa2d' Restart(idx)%axes = 2 diff --git a/ccpp/physics b/ccpp/physics index cf56256643..2ba0850aed 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit cf562566436f5695f318afe5310b3f24e09e9f74 +Subproject commit 2ba0850aed63e48992134f5efe28fcf2a7a545bd diff --git a/ccpp/suites/suite_FV3_GFS_v16_gfdlmpv3.xml b/ccpp/suites/suite_FV3_GFS_v16_gfdlmpv3.xml new file mode 100644 index 0000000000..c02d886b15 --- /dev/null +++ b/ccpp/suites/suite_FV3_GFS_v16_gfdlmpv3.xml @@ -0,0 +1,94 @@ + + + + + + + fv_sat_adj + + + + + GFS_time_vary_pre + GFS_rrtmg_setup + GFS_rad_time_vary + GFS_phys_time_vary + + + + + GFS_rrtmg_pre + GFS_radiation_surface + rad_sw_pre + rrtmg_sw + rrtmg_sw_post + rrtmg_lw + rrtmg_lw_post + GFS_radiation_post + + + + + 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 + satmedmfvdifq + GFS_PBL_generic_post + GFS_GWD_generic_pre + cires_ugwp + cires_ugwp_post + GFS_GWD_generic_post + GFS_suite_stateout_update + + + + + GFS_photochemistry + get_phi_fv3 + GFS_suite_interstitial_3 + GFS_DCNV_generic_pre + samfdeepcnv + GFS_DCNV_generic_post + GFS_SCNV_generic_pre + samfshalcnv + GFS_SCNV_generic_post + GFS_suite_interstitial_4 + cnvc90 + GFS_MP_generic_pre + gfdl_cloud_microphys_v3 + GFS_MP_generic_post + maximum_hourly_diagnostics + GFS_physics_post + + + + + GFS_stochastics + + + + diff --git a/ccpp/suites/suite_FV3_GFS_v17_p8_ugwpv1_tempo.xml b/ccpp/suites/suite_FV3_GFS_v17_p8_ugwpv1_tempo.xml new file mode 100644 index 0000000000..5dfd909e09 --- /dev/null +++ b/ccpp/suites/suite_FV3_GFS_v17_p8_ugwpv1_tempo.xml @@ -0,0 +1,96 @@ + + + + + + + GFS_time_vary_pre + GFS_rrtmg_setup + GFS_rad_time_vary + GFS_phys_time_vary + + + + + GFS_rrtmg_pre + GFS_radiation_surface + rad_sw_pre + rrtmg_sw + rrtmg_sw_post + rrtmg_lw + rrtmg_lw_post + GFS_radiation_post + + + + + 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 + noahmpdrv + sfc_land + sfc_sice + GFS_surface_loop_control_part2 + + + + GFS_surface_composites_post + sfc_diag + sfc_diag_post + GFS_surface_generic_post + GFS_PBL_generic_pre + satmedmfvdifq + GFS_PBL_generic_post + GFS_GWD_generic_pre + ugwpv1_gsldrag + ugwpv1_gsldrag_post + GFS_GWD_generic_post + GFS_suite_stateout_update + + + + + GFS_photochemistry + get_phi_fv3 + GFS_suite_interstitial_3 + GFS_DCNV_generic_pre + samfdeepcnv + GFS_DCNV_generic_post + GFS_SCNV_generic_pre + samfshalcnv + GFS_SCNV_generic_post + GFS_suite_interstitial_4 + cnvc90 + GFS_MP_generic_pre + mp_tempo_pre + + + mp_tempo + + + mp_tempo_post + GFS_MP_generic_post + maximum_hourly_diagnostics + + + + + GFS_stochastics + GFS_physics_post + + + + diff --git a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmpv3_tedmf.xml b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmpv3_tedmf.xml new file mode 100644 index 0000000000..9556c27f9f --- /dev/null +++ b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmpv3_tedmf.xml @@ -0,0 +1,94 @@ + + + + + + + fv_sat_adj + + + + + GFS_time_vary_pre + GFS_rrtmg_setup + GFS_rad_time_vary + GFS_phys_time_vary + + + + + GFS_rrtmg_pre + GFS_radiation_surface + rad_sw_pre + rrtmg_sw + rrtmg_sw_post + rrtmg_lw + rrtmg_lw_post + GFS_radiation_post + + + + + 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 + satmedmfvdifq + GFS_PBL_generic_post + GFS_GWD_generic_pre + unified_ugwp + unified_ugwp_post + GFS_GWD_generic_post + GFS_suite_stateout_update + + + + + GFS_photochemistry + get_phi_fv3 + GFS_suite_interstitial_3 + GFS_DCNV_generic_pre + samfdeepcnv + GFS_DCNV_generic_post + GFS_SCNV_generic_pre + samfshalcnv + GFS_SCNV_generic_post + GFS_suite_interstitial_4 + cnvc90 + GFS_MP_generic_pre + gfdl_cloud_microphys_v3 + GFS_MP_generic_post + maximum_hourly_diagnostics + GFS_physics_post + + + + + GFS_stochastics + + + + diff --git a/ccpp/suites/suite_FV3_coupled_lowres.xml b/ccpp/suites/suite_FV3_coupled_lowres.xml index 8b9a78671a..141fc8df07 100644 --- a/ccpp/suites/suite_FV3_coupled_lowres.xml +++ b/ccpp/suites/suite_FV3_coupled_lowres.xml @@ -19,7 +19,7 @@ rrtmg_sw_post rrtmg_lw rrtmg_lw_post - GFS_rrtmg_post + GFS_radiation_post diff --git a/ccpp/suites/suite_FV3_lowres.xml b/ccpp/suites/suite_FV3_lowres.xml index 5edcef6e28..cdae7a1018 100644 --- a/ccpp/suites/suite_FV3_lowres.xml +++ b/ccpp/suites/suite_FV3_lowres.xml @@ -19,7 +19,7 @@ rrtmg_sw_post rrtmg_lw rrtmg_lw_post - GFS_rrtmg_post + GFS_radiation_post diff --git a/ci/CMakeLists.txt b/ci/CMakeLists.txt index 71412b3afb..43659bb48d 100644 --- a/ci/CMakeLists.txt +++ b/ci/CMakeLists.txt @@ -7,6 +7,15 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMakeModules/Modules) +# Switch to ESMF's cmake files. They have been removed from EMC's modules repo +if(DEFINED ENV{ESMFMKFILE}) + get_filename_component(TEMP_DIR $ENV{ESMFMKFILE} DIRECTORY) + get_filename_component(ESMF_ROOT ${TEMP_DIR} DIRECTORY) + list(APPEND CMAKE_MODULE_PATH ${ESMF_ROOT}/cmake) + include_directories(${ESMF_ROOT}/include) + link_directories(${TEMP_DIR}) +endif() + if(${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 9.0.0) message(FATAL_ERROR "GNU Compiler >= 9 is required") diff --git a/cpl/module_cplfields.F90 b/cpl/module_cplfields.F90 index b23b3b8984..d10eeddaa0 100644 --- a/cpl/module_cplfields.F90 +++ b/cpl/module_cplfields.F90 @@ -160,7 +160,7 @@ module module_cplfields FieldInfo("cpl_scalars ", "s")] ! Import Fields ---------------------------------------- - integer, public, parameter :: NimportFields = 67 + integer, public, parameter :: NimportFields = 67 + 3 + 5 !IVAI: add 3 inst_tracer_diag logical, public :: importFieldsValid(NimportFields) type(ESMF_Field), target, public :: importFields(NimportFields) @@ -183,6 +183,17 @@ module module_cplfields FieldInfo("inst_ice_vis_dir_albedo ", "s"), & FieldInfo("wave_z0_roughness_length ", "s"), & FieldInfo("inst_tracer_diag_aod ", "s"), & +!IVAI: import canopy fields from AQM component + FieldInfo("inst_tracer_diag_claie ", "s"), & + FieldInfo("inst_tracer_diag_cfch ", "s"), & + FieldInfo("inst_tracer_diag_cfrt ", "s"), & + FieldInfo("inst_tracer_diag_cclu ", "s"), & + FieldInfo("inst_tracer_diag_cpopu ", "s"), & +!IVAI: import photolysis diagnostics from AQM component + FieldInfo("inst_tracer_diag_coszens ", "s"), & + FieldInfo("inst_tracer_diag_jo3o1d ", "s"), & + FieldInfo("inst_tracer_diag_jno2 ", "s"), & +!IVAI FieldInfo("ocn_current_zonal ", "s"), & FieldInfo("ocn_current_merid ", "s"), & diff --git a/fv3/atmos_cubed_sphere b/fv3/atmos_cubed_sphere index ac7c36b12a..4f1a5ef4b8 160000 --- a/fv3/atmos_cubed_sphere +++ b/fv3/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit ac7c36b12af82a3ea114825f1272ddb21d29085b +Subproject commit 4f1a5ef4b83461a79a62c8d965ebc2b26cf5eba8 diff --git a/fv3/atmos_model.F90 b/fv3/atmos_model.F90 index 17b5e2bbb2..214319cac8 100644 --- a/fv3/atmos_model.F90 +++ b/fv3/atmos_model.F90 @@ -134,6 +134,8 @@ module atmos_model_mod public atmos_model_get_nth_domain_info public addLsmask2grid public setup_exportdata +public set_fhzero_loop, InitTimeFromIAUOffset +public get_atmos_tracer_types !----------------------------------------------------------------------- ! @@ -270,7 +272,8 @@ subroutine update_atmos_radiation_physics (Atmos) call set_atmosphere_pelist() call mpp_clock_begin(getClock) if (GFS_control%do_skeb) call atmosphere_diss_est (GFS_control%skeb_npass) ! do smoothing for SKEB - call atmos_phys_driver_statein (GFS_Control, GFS_Statein, Atm_block, flip_vc) + ! SA-3D-TKE added GFS_Tbd (kyf) + call atmos_phys_driver_statein (GFS_Control, GFS_Statein, GFS_Tbd, Atm_block, flip_vc) call mpp_clock_end(getClock) !--- if dycore only run, set up the dummy physics output state as the input state @@ -752,7 +755,8 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step) endif ! Populate the GFS_Statein container with the prognostic state ! in Atm_block, which contains the initial conditions/restart data. - call atmos_phys_driver_statein (GFS_control, GFS_statein, Atm_block, flip_vc) + ! SA-3D-TKE added GFS_Tbd (kyf) + call atmos_phys_driver_statein (GFS_control, GFS_statein, GFS_Tbd, Atm_block, flip_vc) ! When asked to calculate 3-dim. tendencies, set Stateout variables to ! Statein variables here in order to capture the first call to dycore @@ -782,33 +786,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step) !--- set the initial diagnostic timestamp diag_time = Time call get_time (Atmos%Time - Atmos%Time_init, sec) - !--- Model should restart at the forecast hours that are multiples of fhzero. - !--- WARNING: For special cases that model needs to restart at non-multiple of fhzero - !--- the fields in first output files are not accumulated from the beginning of - !--- the bucket, but the restart time. - if( GFS_Control%fhzero_array(1) > 0. ) then - fhzero_loop: do i=1,size(GFS_Control%fhzero_array) - tmpflag_fhzero= .false. - if( GFS_Control%fhzero_array(i) > 0.) then - if( i == 1 ) then - if( sec <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true. - else if( i > 1 ) then - if( sec > GFS_Control%fhzero_fhour(i-1)*3600. .and. sec <=GFS_Control%fhzero_fhour(i)*3600. ) & - tmpflag_fhzero = .true. - endif - if( tmpflag_fhzero ) then - GFS_Control%fhzero = GFS_Control%fhzero_array(i) - if( GFS_Control%fhzero > 0) then - sec_lastfhzerofh = (int(sec/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600 - else - sec_lastfhzerofh = 0 - endif - endif - endif - enddo fhzero_loop - else - sec_lastfhzerofh = 0 - endif + call set_fhzero_loop(sec, sec_lastfhzerofh) if (mpp_pe() == mpp_root_pe()) print *,'in atmos_model, fhzero=',GFS_Control%fhzero, 'fhour=',sec/3600.,sec_lastfhzerofh/3600 if (mod((sec-sec_lastfhzerofh),int(GFS_Control%fhzero*3600.)) /= 0) then @@ -856,6 +834,50 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step) end subroutine atmos_model_init ! + !> This will set the forecast hour based on the fhzero_array + !> input. This should handle both an input of a full fhzero + !> array and one that uses increments. + !> + !> @param tmpflag_fhzero logical if current timestep is in between output hours + !> @param[inout] sec time since model initialization, in sec + !> @param[inout] sec_lastfhzerofh time since last fhzero time, in sec + !> + !> @author Daniel Sarmiento @date May 16, 2025 +subroutine set_fhzero_loop(sec, sec_lastfhzerofh) + + logical :: tmpflag_fhzero + integer :: i + integer, intent(inout) :: sec, sec_lastfhzerofh + + !--- Model should restart at the forecast hours that are multiples of fhzero. + !--- WARNING: For special cases that model needs to restart at non-multiple of fhzero + !--- the fields in first output files are not accumulated from the beginning of + !--- the bucket, but the restart time. + if( GFS_Control%fhzero_array(1) > 0. ) then + fhzero_loop: do i=1,size(GFS_Control%fhzero_array) + tmpflag_fhzero= .false. + if( GFS_Control%fhzero_array(i) > 0.) then + if( i == 1 ) then + if( sec <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true. + else if( i > 1 ) then + if( sec > GFS_Control%fhzero_fhour(i-1)*3600. .and. sec <=GFS_Control%fhzero_fhour(i)*3600. ) & + tmpflag_fhzero = .true. + endif + if( tmpflag_fhzero ) then + GFS_Control%fhzero = GFS_Control%fhzero_array(i) + if( GFS_Control%fhzero > 0) then + sec_lastfhzerofh = (int(sec/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600 + else + sec_lastfhzerofh = 0 + endif + endif + endif + enddo fhzero_loop + else + sec_lastfhzerofh = 0 + endif + +end subroutine set_fhzero_loop !####################################################################### ! zero) then - if( time_int - Atmos%iau_offset*3600. > zero ) then - time_int = time_int - Atmos%iau_offset*3600. - else if(seconds == Atmos%iau_offset*3600) then - call get_time (Atmos%Time - diag_time_fhzero, isec_fhzero) - time_int = real(isec_fhzero) - if (mpp_pe() == mpp_root_pe()) write(6,*) "---iseczero",isec_fhzero - endif - endif - time_intfull = real(seconds) - if(Atmos%iau_offset > zero) then - if( time_intfull - Atmos%iau_offset*3600. > zero) then - time_intfull = time_intfull - Atmos%iau_offset*3600. - endif - endif + call InitTimeFromIAUOffset(Atmos, time_int, time_intfull, seconds) if (mpp_pe() == mpp_root_pe()) write(6,*) ' gfs diags time since last bucket empty: ',time_int/3600.,'hrs' call atmosphere_nggps_diag(Atmos%Time) call get_time ( Atmos%Time_step, dtatm_temp) @@ -1028,29 +1036,7 @@ subroutine update_atmos_model_state (Atmos, rc) endif !--- find current fhzero - if( GFS_Control%fhzero_array(1) > 0. ) then - fhzero_loop: do i=1,size(GFS_Control%fhzero_array) - tmpflag_fhzero = .false. - if( GFS_Control%fhzero_array(i) > 0.) then - if( i == 1 ) then - if( seconds <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true. - else if( i > 1 ) then - if( seconds > GFS_Control%fhzero_fhour(i-1)*3600. .and. seconds <= GFS_Control%fhzero_fhour(i)*3600. ) & - tmpflag_fhzero = .true. - endif - if( tmpflag_fhzero) then - GFS_Control%fhzero = GFS_Control%fhzero_array(i) - if( GFS_Control%fhzero > 0) then - sec_lastfhzerofh = (int(seconds/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600 - else - sec_lastfhzerofh = 0 - endif - endif - endif - enddo fhzero_loop - else - sec_lastfhzerofh = 0 - endif + call set_fhzero_loop(seconds,sec_lastfhzerofh) if (mpp_pe() == mpp_root_pe()) print *,'in atmos_model update, fhzero=',GFS_Control%fhzero, 'fhour=',seconds/3600.,sec_lastfhzerofh/3600. if (nint(GFS_Control%fhzero) > 0) then @@ -1078,7 +1064,39 @@ subroutine update_atmos_model_state (Atmos, rc) end subroutine update_atmos_model_state ! + !> This will calculate time if an IAU offest has been defined + !> in the model configuration. + !> + !> @param[inout] atmos the main atmos model configurations + !> @param[inout] time_init model initialization time + !> @param[inout] time_intfull model time remaining + !> @param seconds time since model initialization + !> + !> @author Daniel Sarmiento @date May 16, 2025 + subroutine InitTimeFromIAUOffset(Atmos, time_int, time_intfull, seconds) + + type (atmos_data_type), intent(inout) :: Atmos + real(kind=GFS_kind_phys), intent(inout) :: time_int, time_intfull + integer, intent(inout) :: seconds + integer :: isec_fhzero + + if(Atmos%iau_offset > zero) then + if( time_int - Atmos%iau_offset*3600. > zero ) then + time_int = time_int - Atmos%iau_offset*3600. + else if(seconds == Atmos%iau_offset*3600) then + call get_time (Atmos%Time - diag_time_fhzero, isec_fhzero) + time_int = real(isec_fhzero) + if (mpp_pe() == mpp_root_pe()) write(6,*) "---iseczero",isec_fhzero + endif + endif + time_intfull = real(seconds) + if(Atmos%iau_offset > zero) then + if( time_intfull - Atmos%iau_offset*3600. > zero) then + time_intfull = time_intfull - Atmos%iau_offset*3600. + endif + endif + end subroutine InitTimeFromIAUOffset !####################################################################### ! @@ -1340,8 +1358,12 @@ subroutine update_atmos_chemistry(state, rc) real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: q +!IVAI: add coszens, jo3o1d, jno2, claie, cfch, cfrt, cclu, cpopu real(ESMF_KIND_R8), dimension(:,:), pointer :: aod, area, canopy, cmm, & - dqsfc, dtsfc, fice, flake, focn, fsnow, hpbl, nswsfc, oro, psfc, & + claie, cfch, cfrt, cclu, cpopu, & !IVAI + dqsfc, dtsfc, fice, flake, focn, fsnow, hpbl, & + coszens, jo3o1d, jno2, & !IVAI + nswsfc, oro, psfc, & q2m, rain, rainc, rca, shfsfc, slmsk, stype, swet, t2m, tsfc, & u10m, uustar, v10m, vfrac, xlai, zorl, vtype @@ -1368,6 +1390,44 @@ subroutine update_atmos_chemistry(state, rc) call cplFieldGet(state,'inst_tracer_diag_aod', farrayPtr2d=aod, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__, rcToReturn=rc)) return + +!IVAI: case ('import') canopy arrays read in via 'aqm_emis_read' + + if (GFS_control%do_canopy) then + call cplFieldGet(state,'inst_tracer_diag_claie', farrayPtr2d=claie, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_cfch', farrayPtr2d=cfch, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_cfrt', farrayPtr2d=cfrt, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_cclu', farrayPtr2d=cclu, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_cpopu', farrayPtr2d=cpopu, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + end if + +!IVAI: case ('import') photdiag arrays + call cplFieldGet(state,'inst_tracer_diag_coszens', farrayPtr2d=coszens, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_jo3o1d', farrayPtr2d=jo3o1d, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call cplFieldGet(state,'inst_tracer_diag_jno2', farrayPtr2d=jno2, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return +!IVAI end if !--- do not import tracer concentrations by default @@ -1442,6 +1502,123 @@ subroutine update_atmos_chemistry(state, rc) GFS_IntDiag%aod(im) = aod(i,j) enddo enddo + + if (GFS_control%do_canopy) then +!IVAI: case ('import') canopy arrays read in via aqm_emis_read +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, claie) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%claie(im) = claie(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, cfch) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%cfch(im) = cfch(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, cfrt) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%cfrt(im) = cfrt(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, cclu) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%cclu(im) = cclu(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, cpopu) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%cpopu(im) = cpopu(i,j) + enddo + enddo + endif ! GFS_control%do_canopy + +!IVAI: case ('import') photdiag arrays +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, coszens) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%coszens(im) = coszens(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, jo3o1d) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%jo3o1d(im) = jo3o1d(i,j) + enddo + enddo + +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, GFS_Control, GFS_Intdiag, jno2) & +!$OMP private (j, jb, i, ib, nb, ix, im) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + im = GFS_Control%chunk_begin(nb)+ix-1 + GFS_IntDiag%jno2(im) = jno2(i,j) + enddo + enddo +!IVAI end if if (GFS_control%debug) then @@ -1450,6 +1627,33 @@ subroutine update_atmos_chemistry(state, rc) if (GFS_control%cplaqm) & write(6,'("update_atmos: ",a,": aod - min/max ",3g16.6)') & trim(state), minval(aod), maxval(aod) +!IVAI: case ('import') canopy arrays read via aqm_emis_read + if (GFS_control%cplaqm .and. GFS_control%do_canopy) & + write(6,'("update_atmos: ",a,": claie - min/max ",3g16.6)') & + trim(state), minval(claie), maxval(claie) + if (GFS_control%cplaqm .and. GFS_control%do_canopy) & + write(6,'("update_atmos: ",a,": cfch - min/max ",3g16.6)') & + trim(state), minval(cfch), maxval(cfch) + if (GFS_control%cplaqm .and. GFS_control%do_canopy) & + write(6,'("update_atmos: ",a,": cfrt - min/max ",3g16.6)') & + trim(state), minval(cfrt), maxval(cfrt) + if (GFS_control%cplaqm .and. GFS_control%do_canopy) & + write(6,'("update_atmos: ",a,": cclu - min/max ",3g16.6)') & + trim(state), minval(cclu), maxval(cclu) + if (GFS_control%cplaqm .and. GFS_control%do_canopy) & + write(6,'("update_atmos: ",a,": cpopu - min/max ",3g16.6)') & + trim(state), minval(cpopu), maxval(cpopu) +!IVAI: case ('import') photdiag arrays + if (GFS_control%cplaqm) & + write(6,'("update_atmos: ",a,": coszens - min/max ",3g16.6)') & + trim(state), minval(coszens), maxval(coszens) + if (GFS_control%cplaqm) & + write(6,'("update_atmos: ",a,": jo3o1d - min/max ",3g16.6)') & + trim(state), minval(jo3o1d), maxval(jo3o1d) + if (GFS_control%cplaqm) & + write(6,'("update_atmos: ",a,": jno2 - min/max ",3g16.6)') & + trim(state), minval(jno2), maxval(jno2) +!IVAI end if case ('export') diff --git a/fv3/fv3_cap.F90 b/fv3/fv3_cap.F90 index 49f13d0d90..56a0c2ba46 100644 --- a/fv3/fv3_cap.F90 +++ b/fv3/fv3_cap.F90 @@ -9,7 +9,8 @@ ! 18 Apr 2017: J. Wang set up fcst grid component and write grid components ! 24 Jul 2017: J. Wang initialization and time stepping changes for coupling ! 02 Nov 2017: J. Wang Use Gerhard's transferable RouteHandle -! +! 20 May 2025: D. Sarmiento Handle output hour array in seperate subroutines +! module fv3atm_cap_mod @@ -53,6 +54,7 @@ module fv3atm_cap_mod implicit none private public SetServices + public OutputHours_FrequencyInput, OutputHours_ArrayInput ! !----------------------------------------------------------------------- ! @@ -943,27 +945,9 @@ subroutine InitializeAdvertise(gcomp, rc) call ESMF_ConfigGetAttribute(CF,valueList=outputfh2,label='output_fh:', & count=noutput_fh, rc=rc) if(outputfh2(2) == -1) then - !--- output_hf is output frequency, the second item is -1 + !--- output_fh is output frequency, the second item is -1 lfreq = .true. - nfh = 0 - if( nfhmax>output_startfh) nfh = nint((nfhmax-output_startfh)/outputfh2(1)) + 1 - if( nfh > 0) then - allocate(output_fh(nfh)) - if( output_startfh == 0) then - output_fh(1) = dt_atmos/3600. - else - output_fh(1) = output_startfh - endif - do i=2,nfh - output_fh(i) = (i-1)*outputfh2(1) + output_startfh - ! Except fh000, which is the first time output, if any other of the - ! output time is not integer hour, set lflname_fulltime to be true, so the - ! history file names will contain the full time stamp (HHH-MM-SS). - if(.not.lflname_fulltime) then - if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. - endif - enddo - endif + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) endif endif if( noutput_fh /= 2 .or. .not. lfreq ) then @@ -972,32 +956,7 @@ subroutine InitializeAdvertise(gcomp, rc) call ESMF_ConfigGetAttribute(CF,valueList=output_fh,label='output_fh:', & count=noutput_fh, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - if( output_startfh == 0) then - ! If the output time in output_fh array contains first time stamp output, - ! check the rest of output time, otherwise, check all the output time. - ! If any of them is not integer hour, the history file names will - ! contain the full time stamp (HHH-MM-SS) - ist = 1 - if(output_fh(1)==0) then - output_fh(1) = dt_atmos/3600. - ist= 2 - endif - do i=ist,noutput_fh - if(.not.lflname_fulltime) then - if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. - endif - enddo - else - do i=1,noutput_fh - output_fh(i) = output_startfh + output_fh(i) - ! When output_startfh >0, check all the output time, if any of - ! them is not integer hour, set lflname_fulltime to be true. The - ! history file names will contain the full time stamp (HHH-MM-SS). - if(.not.lflname_fulltime) then - if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. - endif - enddo - endif + call OutputHours_ArrayInput(noutput_fh,output_startfh) endif endif ! end loutput_fh endif @@ -1030,6 +989,80 @@ subroutine InitializeAdvertise(gcomp, rc) end subroutine InitializeAdvertise !----------------------------------------------------------------------------- + !> This will calculate output hours if the user has stated a + !> an fhzero frequency. + !> + !> @param[inout] nfhmax maximum number of forecast hours + !> @param[inout] output_startfh ouptut start time + !> @param[inout] outputfh2 user defined forecast hour configuration + !> + !> @author Daniel Sarmiento @date May 16, 2025 + subroutine OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + integer :: nfh, i + real, intent(inout) :: nfhmax, output_startfh, outputfh2(2) + + nfh = 0 + if( nfhmax>output_startfh) nfh = nint((nfhmax-output_startfh)/outputfh2(1)) + 1 + if( nfh > 0) then + allocate(output_fh(nfh)) + if( output_startfh == 0) then + output_fh(1) = dt_atmos/3600. + else + output_fh(1) = output_startfh + endif + do i=2,nfh + output_fh(i) = (i-1)*outputfh2(1) + output_startfh + ! Except fh000, which is the first time output, if any other of the + ! output time is not integer hour, set lflname_fulltime to be true, so the + ! history file names will contain the full time stamp (HHH-MM-SS). + if(.not.lflname_fulltime) then + if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. + endif + enddo + endif + end subroutine OutputHours_FrequencyInput + + !> This will calculate output hours if the user has stated a + !> an array of desired output hours. + !> + !> @param[inout] noutput_fh index of output hours array + !> @param[inout] output_startfh ouptut start time + !> + !> @author Daniel Sarmiento @date May 16, 2025 + subroutine OutputHours_ArrayInput(noutput_fh,output_startfh) + + integer :: ist, i + integer, intent(inout) :: noutput_fh + real, intent(inout) :: output_startfh + + if( output_startfh == 0) then + ! If the output time in output_fh array contains first time stamp output, + ! check the rest of output time, otherwise, check all the output time. + ! If any of them is not integer hour, the history file names will + ! contain the full time stamp (HHH-MM-SS) + ist = 1 + if(output_fh(1)==0) then + output_fh(1) = dt_atmos/3600. + ist= 2 + endif + do i=ist,noutput_fh + if(.not.lflname_fulltime) then + if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. + endif + enddo + else + do i=1,noutput_fh + output_fh(i) = output_startfh + output_fh(i) + ! When output_startfh >0, check all the output time, if any of + ! them is not integer hour, set lflname_fulltime to be true. The + ! history file names will contain the full time stamp (HHH-MM-SS). + if(.not.lflname_fulltime) then + if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true. + endif + enddo + endif + + end subroutine OutputHours_ArrayInput subroutine InitializeRealize(gcomp, rc) type(ESMF_GridComp) :: gcomp diff --git a/fv3/io/module_wrt_grid_comp.F90 b/fv3/io/module_wrt_grid_comp.F90 index 7b2c1dd499..9baea58095 100644 --- a/fv3/io/module_wrt_grid_comp.F90 +++ b/fv3/io/module_wrt_grid_comp.F90 @@ -59,6 +59,7 @@ module module_wrt_grid_comp !----------------------------------------------------------------------- private + public get_outfile, lambert, rtll ! !----------------------------------------------------------------------- ! diff --git a/fv3/module_fcst_grid_comp.F90 b/fv3/module_fcst_grid_comp.F90 index aa8241ffe8..cea2dfdf2c 100644 --- a/fv3/module_fcst_grid_comp.F90 +++ b/fv3/module_fcst_grid_comp.F90 @@ -759,57 +759,10 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) ! Time_step = set_time (dt_atmos,0) if (mype == 0) write(*,*)'time_init=', date_init,'time=',date,'time_end=',date_end,'dt_atmos=',dt_atmos - -! set up forecast time array that controls when to write out restart files - frestart = 0 - call get_time(Time_end - Time_init, total_inttime) -! set iau offset time - Atmos%iau_offset = iau_offset - if(iau_offset > 0 ) then - iautime = set_time(iau_offset * 3600, 0) - endif -! if the second item is -1, the first number is frequency - freq_restart = .false. - if(num_restart_interval == 2) then - if(restart_interval(2)== -1) freq_restart = .true. - endif - if(freq_restart) then - if(restart_interval(1) >= 0) then - tmpvar = restart_interval(1) * 3600 - Time_step_restart = set_time (tmpvar, 0) - if(iau_offset > 0 ) then - Time_restart = Time_init + iautime + Time_step_restart - frestart(1) = tmpvar + iau_offset *3600 - else - Time_restart = Time_init + Time_step_restart - frestart(1) = tmpvar - endif - if(restart_interval(1) > 0) then - i = 2 - do while ( Time_restart < Time_end ) - frestart(i) = frestart(i-1) + tmpvar - Time_restart = Time_restart + Time_step_restart - i = i + 1 - enddo - endif - endif -! otherwise it is an array with forecast time at which the restart files will be written out - else if(num_restart_interval >= 1) then - if(num_restart_interval == 1 .and. restart_interval(1) == 0 ) then - frestart(1) = total_inttime - else - if(iau_offset > 0 ) then - restart_starttime = iau_offset *3600 - else - restart_starttime = 0 - endif - do i=1,num_restart_interval - frestart(i) = restart_interval(i) * 3600. + restart_starttime - enddo - endif - endif -! if to write out restart at the end of forecast - if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'total_inttime=',total_inttime + + call fcst_time_array_setup(Time_init, Time_end, Time_step_restart, & + Time_restart, num_restart_interval, & + restart_interval) !------ initialize component models ------ @@ -1197,6 +1150,83 @@ subroutine create_bundle_and_add_it_to_state(name_fb, state, rc) end subroutine create_bundle_and_add_it_to_state end subroutine fcst_initialize + + !> Create forecast hour time array. This will be used + !> to dictate when restart files are going to be written. + !> + !> @param[inout] Time_init model initialization time + !> @param[inout] Time_end model end time + !> @param[inout] Time_step_restart restart time based on restart_interval + !> @param[inout] Time_restart calculated restart time + !> @param[inout] num_restart_interval user defined restart interval + !> @param[inout] restart_interval restart interval, allocatable + !> + !> @author Daniel Sarmiento @date May 16, 2025 + subroutine fcst_time_array_setup(Time_init, Time_end, Time_step_restart, & + Time_restart, num_restart_interval, & + restart_interval) + + type(time_type), intent(inout) :: Time_init, Time_end, & + Time_step_restart, & + Time_restart + type(time_type) :: iautime + integer, intent(inout) :: num_restart_interval + integer :: total_inttime, tmpvar, & + i, restart_starttime + logical :: freq_restart + real, dimension(:), allocatable, intent(inout) :: restart_interval + + ! set up forecast time array that controls when to write out restart files + frestart = 0 + call get_time(Time_end - Time_init, total_inttime) + ! set iau offset time + Atmos%iau_offset = iau_offset + if(iau_offset > 0 ) then + iautime = set_time(iau_offset * 3600, 0) + endif + ! if the second item is -1, the first number is frequency + freq_restart = .false. + if(num_restart_interval == 2) then + if(restart_interval(2)== -1) freq_restart = .true. + endif + if(freq_restart) then + if(restart_interval(1) >= 0) then + tmpvar = restart_interval(1) * 3600 + Time_step_restart = set_time (tmpvar, 0) + if(iau_offset > 0 ) then + Time_restart = Time_init + iautime + Time_step_restart + frestart(1) = tmpvar + iau_offset *3600 + else + Time_restart = Time_init + Time_step_restart + frestart(1) = tmpvar + endif + if(restart_interval(1) > 0) then + i = 2 + do while ( Time_restart < Time_end ) + frestart(i) = frestart(i-1) + tmpvar + Time_restart = Time_restart + Time_step_restart + i = i + 1 + enddo + endif + endif + ! otherwise it is an array with forecast time at which the restart files will be written out + else if(num_restart_interval >= 1) then + if(num_restart_interval == 1 .and. restart_interval(1) == 0 ) then + frestart(1) = total_inttime + else + if(iau_offset > 0 ) then + restart_starttime = iau_offset *3600 + else + restart_starttime = 0 + endif + do i=1,num_restart_interval + frestart(i) = restart_interval(i) * 3600. + restart_starttime + enddo + endif + endif + ! if to write out restart at the end of forecast + if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'total_inttime=',total_inttime + end subroutine fcst_time_array_setup ! !----------------------------------------------------------------------- !####################################################################### diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5667d3d619..1d00b4b764 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,11 @@ # This file sets up CTest unit tests for fv3atm. # # Alex Richert, Jan 2024 +# +# Add More unit tests +# Daniel Sarmiento, Jun 2025 + +find_package(PIO REQUIRED COMPONENTS C Fortran) # Stage test data execute_process(COMMAND cmake -E create_symlink @@ -10,19 +15,19 @@ execute_process(COMMAND cmake -E create_symlink function(add_fv3atm_mpi_test TESTNAME) add_executable(${TESTNAME} ${TESTNAME}.F90) - target_link_libraries(${TESTNAME} PRIVATE fv3atm MPI::MPI_Fortran PIO::PIO_Fortran) + target_link_libraries(${TESTNAME} PRIVATE ufsatm_fv3 MPI::MPI_Fortran PIO::PIO_C) add_test(${TESTNAME} ${MPIEXEC_EXECUTABLE} -n 2 ${CMAKE_CURRENT_BINARY_DIR}/${TESTNAME}) endfunction() function(add_fv3atm_serial_test TESTNAME) add_executable(${TESTNAME} ${TESTNAME}.F90) - target_link_libraries(${TESTNAME} PRIVATE fv3atm) + target_link_libraries(${TESTNAME} PRIVATE ufsatm_fv3) add_test(${TESTNAME} ${CMAKE_CURRENT_BINARY_DIR}/${TESTNAME}) endfunction() -#foreach(testname test_write_netcdf) -# add_fv3atm_mpi_test(${testname}) -#endforeach() +foreach(testname test_atmos_model test_fv3_cap test_module_wrt_grid_comp) + add_fv3atm_mpi_test(${testname}) +endforeach() foreach(testname test_post_nems_routines) add_fv3atm_serial_test(${testname}) diff --git a/tests/test_atmos_model.F90 b/tests/test_atmos_model.F90 new file mode 100644 index 0000000000..1b79838bae --- /dev/null +++ b/tests/test_atmos_model.F90 @@ -0,0 +1,192 @@ +program test_atmos_model + use atmos_model_mod, only: set_fhzero_loop, InitTimeFromIAUOffset, & + get_atmos_tracer_types, atmos_data_type + use GFS_typedefs, only: GFS_control_type, GFS_kind_phys => kind_phys + use time_manager_mod, only: time_type, set_time, get_time, operator(-) + use tracer_manager_mod, only: get_number_tracers + use field_manager_mod, only: MODEL_ATMOS + use mpp_mod, only: mpp_init, mpp_exit, FATAL, mpp_error + use CCPP_data, only: GFS_control + + implicit none + + integer :: test_passed, total_tests + integer :: suite_passed, suite_total + + ! Test Suite 1: set_fhzero_loop + call test_set_fhzero_loop_suite() + + ! Test Suite 2: get_atmos_tracer_types + call test_get_atmos_tracer_types_suite() + +contains + + !============================================================================ + ! TEST SUITE 1: set_fhzero_loop + !============================================================================ + subroutine test_set_fhzero_loop_suite() + integer :: sec, sec_lastfhzerofh + + ! Test 1: Basic functionality with single fhzero value + call test_single_fhzero() + + ! Test 2: Multiple fhzero array values + call test_multiple_fhzero() + + ! Test 3: Edge case with zero and negative values + call test_fhzero_edge_cases() + + end subroutine test_set_fhzero_loop_suite + + subroutine test_single_fhzero() + integer :: sec, sec_lastfhzerofh + + ! Setup + GFS_control%fhzero_array(1) = 6.0_GFS_kind_phys + GFS_control%fhzero_fhour(1) = 24.0_GFS_kind_phys + + ! Test case: Time within first interval + sec = 10800 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 6.0_GFS_kind_phys .or. sec_lastfhzerofh /= 0) then + print *, "Incorrect handling of single fhzero value" + stop 1 + end if + + end subroutine test_single_fhzero + + subroutine test_multiple_fhzero() + integer :: sec, sec_lastfhzerofh + + ! Setup + GFS_control%fhzero_array = [3.0_GFS_kind_phys, 6.0_GFS_kind_phys] + GFS_control%fhzero_fhour = [12.0_GFS_kind_phys, 24.0_GFS_kind_phys] + + ! Test first interval + sec = 7200 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 3.0_GFS_kind_phys) then + print *, "Incorrect handling of fhzero array" + stop 2 + end if + + ! Test second interval + sec = 50400 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 6.0_GFS_kind_phys) then + print *, "Incorrect handling of fhzero array" + stop 3 + end if + + end subroutine test_multiple_fhzero + + subroutine test_fhzero_edge_cases() + integer :: sec, sec_lastfhzerofh + + ! Test zero fhzero value + GFS_control%fhzero_array = [0.0_GFS_kind_phys, 6.0_GFS_kind_phys] + GFS_control%fhzero_fhour = [6.0_GFS_kind_phys, 12.0_GFS_kind_phys] + + sec = 3600 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (sec_lastfhzerofh /= 0) then + print *, "Incorrect handling of fh = 0 case" + stop 4 + end if + + end subroutine test_fhzero_edge_cases + + !============================================================================ + ! TEST SUITE 2: get_atmos_tracer_types + !============================================================================ + subroutine test_get_atmos_tracer_types_suite() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Test 1: Basic functionality with mock tracers + call test_tracer_basic_functionality() + + ! Test 2: Test with chemistry tracers + call test_chemistry_tracers() + + ! Test 3: Edge cases + call test_tracer_edge_cases() + + end subroutine test_get_atmos_tracer_types_suite + + subroutine test_tracer_basic_functionality() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! For this test, we'll simulate having 5 tracers + num_tracers = 5 + allocate(tracer_types(num_tracers)) + + ! Initialize all to zero (default) + tracer_types = 0 + + if (any(tracer_types /= 0)) then + print *, "Tracer type array being rewritten" + stop 5 + end if + + deallocate(tracer_types) + end subroutine test_tracer_basic_functionality + + subroutine test_chemistry_tracers() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Simulate having tracers with chemistry types + num_tracers = 8 + allocate(tracer_types(num_tracers)) + + ! Manually set tracer types to simulate: + ! Tracers 1-3: generic (0) + ! Tracers 4-6: chemistry prognostic (1) + ! Tracers 7-8: chemistry diagnostic (2) + tracer_types = [0, 0, 0, 1, 1, 1, 2, 2] + + ! Test generic tracers are contiguous + if (any(tracer_types(1:3) /= 0)) then + print *, "Tracer type array being rewritten or rearranged" + stop 6 + end if + + ! Test prognostic tracers are contiguous + if (any(tracer_types(4:6) /= 1)) then + print *, "Tracer type array being rewritten or rearranged" + stop 7 + end if + + ! Test diagnostic tracers are contiguous + if (any(tracer_types(7:8) /= 2)) then + print *, "Tracer type array being rewritten or rearranged" + stop 8 + end if + + deallocate(tracer_types) + end subroutine test_chemistry_tracers + + subroutine test_tracer_edge_cases() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Test with large number of tracers + num_tracers = 100 + allocate(tracer_types(num_tracers)) + tracer_types = 0 + + if (size(tracer_types) /= 100) then + print *, "Tracer type array missing values when array is large" + stop 9 + end if + + deallocate(tracer_types) + end subroutine test_tracer_edge_cases + +end program test_atmos_model diff --git a/tests/test_fv3_cap.F90 b/tests/test_fv3_cap.F90 new file mode 100644 index 0000000000..7d98fb5717 --- /dev/null +++ b/tests/test_fv3_cap.F90 @@ -0,0 +1,229 @@ +program test_output_hours + use fv3atm_cap_mod, only: OutputHours_FrequencyInput, OutputHours_ArrayInput + use module_fv3_config, only: dt_atmos, output_fh + use module_fv3_io_def, only: lflname_fulltime + + implicit none + + ! Test variables + integer :: test_passed, test_failed + integer :: i, expected_size + logical :: test_result + + ! Variables for testing + real :: nfhmax, output_startfh, outputfh2(2) + integer :: noutput_fh + + dt_atmos = 1800 + + call test_frequency_input() + call test_array_input() + +contains + + ! Test OutputHours_FrequencyInput subroutine + subroutine test_frequency_input() + + !============================================ + ! Test 1: Basic frequency input with start at 0 + nfhmax = 24.0 + output_startfh = 0.0 + outputfh2(1) = 3.0 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + expected_size = 9 + + if (size(output_fh) /= expected_size) then + !Expected size should be 9 + print *, "Size of generated output_fh is incorrect" + stop 1 + end if + if (abs(output_fh(1) - dt_atmos/3600.0) > 1e-6) then + !First output should be dt_atmos/3600 + print *, "First output time is incorrect" + stop 2 + end if + ! lflname_fulltime should be false + if (lflname_fulltime) then + !Excluding first element, output_fh are all integers + print *, "lflname_fulltime bool set incorrectly" + stop 3 + end if + + !============================================ + ! Test 2: Frequency input with non-zero start + nfhmax = 48.0 + output_startfh = 6.0 + outputfh2(1) = 6.0 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + expected_size = 8 + + !Expected size should be 8 + if (size(output_fh) /= expected_size) then + print *, "Size of generated output_fh is incorrect" + stop 4 + end if + + ! First value (should be 6.0) + if (abs(output_fh(1) - 6.0) > 1e-6) then + print *, "First output time is incorrect" + stop 5 + end if + + ! lflname_fulltime should be false + if (lflname_fulltime) then + !Excluding first element, output_fh are all integers + print *, "lflname_fulltime bool set incorrectly" + stop 6 + end if + + !============================================ + ! Test 3: Frequency that creates non-integer hours + ! Only checking lflname_fulltime since other aspects of + ! array generation were already checked. + nfhmax = 10.0 + output_startfh = 0.0 + outputfh2(1) = 2.5 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + ! Check lflname_fulltime, should be True since non-integer values exist + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 7 + end if + + !============================================ + ! Test 4: nfhmax equals output_startfh + nfhmax = 6.0 + output_startfh = 6.0 + outputfh2(1) = 3.0 + outputfh2(2) = -1.0 + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + ! output_fh should not allocate when nfhmax == output_startfh + if (allocated(output_fh)) then + print *, "output_fh was allocated when output start time was equal to nfmax" + stop 8 + end if + + end subroutine test_frequency_input + + subroutine test_array_input() + + !============================================ + ! Test 1: Basic array input with start at 0 + noutput_fh = 5 + output_startfh = 0.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 3.0, 6.0, 9.0, 12.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check first value. Should be dt_atmos/3600 + if (abs(output_fh(1) - dt_atmos/3600.0) > 1e-6) then + print *, "First output time is incorrect" + stop 9 + end if + + ! Check lflname_fulltime, should be false + !Excluding first element, output_fh are all integers + if (lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 10 + end if + + !============================================ + ! Test 2: Array input with non-zero start + noutput_fh = 4 + output_startfh = 6.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 6.0, 12.0, 18.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check values (should be shifted by 6: 6, 12, 18, 24) + if (abs(output_fh(1) - 6.0) > 1e-6 .or. & + abs(output_fh(2) - 12.0) > 1e-6 .or. & + abs(output_fh(3) - 18.0) > 1e-6 .or. & + abs(output_fh(4) - 24.0) > 1e-6) then + print *, "output_fh array values were not shifted correctly or not allocated correctly" + stop 11 + end if + + ! Check lflname_fulltime (should be false) + if (lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 12 + end if + + !============================================ + ! Test 3: Array with non-integer hours + noutput_fh = 4 + output_startfh = 0.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 1.5, 3.0, 4.5, 6.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + test_result = .true. + print *, ' Output hours:', output_fh + + ! Check lflname_fulltime (should be true) + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 13 + end if + + !============================================ + ! Test 4: Array with fractional hours from non-zero start + noutput_fh = 3 + output_startfh = 0.5 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 3.0, 6.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check values (should be 0.5, 3.5, 6.5) + if (abs(output_fh(1) - 0.5) > 1e-6 .or. & + abs(output_fh(2) - 3.5) > 1e-6 .or. & + abs(output_fh(3) - 6.5) > 1e-6) then + print *, "output_fh array values were not shifted correctly or not allocated correctly" + stop 14 + end if + + ! Check lflname_fulltime (should be true) + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 15 + end if + + end subroutine test_array_input + +end program test_output_hours diff --git a/tests/test_module_wrt_grid_comp.F90 b/tests/test_module_wrt_grid_comp.F90 new file mode 100644 index 0000000000..a3c67abad7 --- /dev/null +++ b/tests/test_module_wrt_grid_comp.F90 @@ -0,0 +1,181 @@ +program test_module_wrt_grid_comp + use module_wrt_grid_comp, only: get_outfile, lambert, rtll + implicit none + + call test_get_outfile() + call test_lambert() + call test_rtll() + +contains + + !--------------------------------------------------------------------------- + ! Test get_outfile subroutine + !--------------------------------------------------------------------------- + subroutine test_get_outfile() + character(len=128) :: filename(2000,3) + character(len=128) :: outfile_name(100) + integer :: noutfile + character(len=100) :: test_name + + ! Test 1: Basic test with unique filenames + filename = '' + filename(1,1) = 'file1.nc' + filename(2,1) = 'file2.nc' + filename(1,2) = 'file3.nc' + filename(1,3) = 'file4.nc' + + call get_outfile(3, filename, outfile_name, noutfile) + if ( trim(outfile_name(1)) /= "file1.nc" .or. & + trim(outfile_name(2)) /= "file2.nc" .or. & + trim(outfile_name(3)) /= "file3.nc" .or. & + trim(outfile_name(4)) /= "file4.nc") then + print *, "Filename trim function not working as intended" + stop 1 + + end if + + end subroutine test_get_outfile + + !--------------------------------------------------------------------------- + ! Test lambert subroutine + !--------------------------------------------------------------------------- + subroutine test_lambert() + real(8) :: stlat1, stlat2, c_lat, c_lon + real(8) :: glon, glat, x, y + real(8) :: glon_inv, glat_inv, x_out, y_out + real(8) :: true_x, true_y + real(8), parameter :: tol = 1.0e2 ! Difference tolerance + + ! Test 1: Forward transformation (glon,glat) -> (x,y) + ! National mall, Washington DC + stlat1 = 38.5_8 + stlat2 = 39.5_8 + c_lat = 39.0_8 + c_lon = -77.0_8 + glon = -77.0353_8 + glat = 38.8895_8 + + true_x = -3055.18 + true_y = -12286.37 + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( abs(true_x - x) > tol .or. abs(true_y - y) > tol ) then + print *, "Forward transformation incorrect" + stop 2 + end if + + ! Test inverse transformation + glon_inv = 0.0_8 + glat_inv = 0.0_8 + call lambert(stlat1, stlat2, c_lat, c_lon, glon_inv, glat_inv, x, y, -1) + + if ( (glon - glon_inv) > tol .or. (glat - glat_inv) > tol ) then + print *, "Results from glon,glat -> lon,lat -> glon,glat not identical" + stop 3 + end if + + ! Test 2: Special case where stlat1 == stlat2 + ! Italy + stlat1 = 45.5_8 + stlat2 = 45.5_8 + c_lat = 45.5_8 + c_lon = 11.0_8 + glat = 45.4642_8 + glon = 9.1900_8 + + true_x = -141149.15 + true_y = -2390.66 + + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( (true_x - x) > tol .or. (true_y - y) > tol ) then + print *, "Issue when stlat1 and stlat2 are equal" + stop 4 + end if + + ! Test 3: Point at projection center + stlat1 = 38.5_8 + stlat2 = 39.5_8 + c_lat = 39.0_8 + c_lon = -77.0_8 + glon = -77.0_8 + glat = 39.0_8 + + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( abs(x) > 1e-3 .or. abs(y) > 1e-3 ) then + print *, "Issue when point is at the projection center" + stop 5 + end if + + end subroutine test_lambert + + !--------------------------------------------------------------------------- + ! Test rtll subroutine + !--------------------------------------------------------------------------- + subroutine test_rtll() + real(8) :: tlmd, tphd, almd, aphd, tlm0d, tph0d + real(8) :: true_almd, true_aphd + real(8), parameter :: tol = 1.0e-0 + + ! Test 1: No rotation + tlm0d = 0.0_8 + tph0d = 0.0_8 + tlmd = 0.0_8 + tphd = 0.0_8 + true_almd = 0.0_8 + true_aphd = 0.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Coordinates were rotated when no rotation was expected" + stop 6 + end if + + ! Test 2: North pole rotation + tlm0d = 0.0_8 + tph0d = 90.0_8 + tlmd = 0.0_8 + tphd = 0.0_8 + true_almd = 0.0_8 + true_aphd = 90.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Incorrect rotation when 90 deg rotation was expected" + stop 7 + end if + + ! Test 3: Central Europe + tlm0d = -170.0_8 + tph0d = 40.0_8 + tlmd = 10.0_8 + tphd = 50.0_8 + true_almd = -76.16_8 + true_aphd = 83.57_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Incorrect rotation of coordiantes" + stop 8 + end if + + ! Test 4: Longitude wrapping + tlm0d = 0.0_8 + tph0d = 45.0_8 + tlmd = 179.0_8 + tphd = 0.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( almd > 180 .or. almd < -180) then + print *, "Incorrect rotation when longitude is at the edge of the domain" + stop 9 + end if + + end subroutine test_rtll + +end program test_module_wrt_grid_comp