diff --git a/CMakeLists.txt b/CMakeLists.txt index 125c17fa9..db576a054 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -86,7 +86,9 @@ add_library(fv3atm cpl/module_cplfields.F90 cpl/module_cap_cpl.F90 io/FV3GFS_io.F90 + io/FV3GFS_restart_io.F90 io/module_write_netcdf.F90 + io/module_write_restart_netcdf.F90 io/module_fv3_io_def.F90 io/module_write_internal_state.F90 io/module_wrt_grid_comp.F90 diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere index a4450fbed..9d5bed8e9 160000 --- a/atmos_cubed_sphere +++ b/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit a4450fbed77aca481c50db1623454fe7ff58d8bd +Subproject commit 9d5bed8e932a188d714c8c8a770e44801dab4750 diff --git a/atmos_model.F90 b/atmos_model.F90 index 2a5091bd6..e6bbb72c7 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -54,7 +54,7 @@ module atmos_model_mod use fms_mod, only: check_nml_error use diag_manager_mod, only: diag_send_complete_instant use time_manager_mod, only: time_type, get_time, get_date, & - operator(+), operator(-),real_to_time_type + operator(+), operator(-), real_to_time_type use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_number_tracers, get_tracer_names, & get_tracer_index, NO_TRACER @@ -94,9 +94,15 @@ module atmos_model_mod FV3GFS_GFS_checksum, & FV3GFS_diag_register, FV3GFS_diag_output, & DIAG_SIZE +use FV3GFS_restart_io_mod, only: FV3GFS_restart_register, & + fv_phy_restart_output, & + fv_sfc_restart_output +use fv_ufs_restart_io_mod, only: fv_dyn_restart_register, & + fv_dyn_restart_output use fv_iau_mod, only: iau_external_data_type,getiauforcing,iau_initialize use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout, & - output_fh, fcst_mpi_comm, fcst_ntasks + output_fh, fcst_mpi_comm, fcst_ntasks, & + quilting_restart use module_block_data, only: block_atmos_copy, block_data_copy, & block_data_copy_or_fill, & block_data_combine_fractions @@ -177,7 +183,6 @@ module atmos_model_mod logical :: debug = .false. !logical :: debug = .true. logical :: sync = .false. -logical :: restart_endfcst = .false. real :: avg_max_length=3600. logical :: ignore_rst_cksum = .false. namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, ccpp_suite, avg_max_length, & @@ -735,6 +740,10 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step) call GFS_restart_populate (GFS_restart_var, GFS_control, GFS_data%Statein, GFS_data%Stateout, GFS_data%Sfcprop, & GFS_data%Coupling, GFS_data%Grid, GFS_data%Tbd, GFS_data%Cldprop, GFS_data%Radtend, & GFS_data%IntDiag, Init_parm, GFS_Diag) + if (quilting_restart) then + call fv_dyn_restart_register (Atm(mygrid)) + call FV3GFS_restart_register (GFS_data%Sfcprop, GFS_restart_var, Atm_block, GFS_control) + endif call FV3GFS_restart_read (GFS_data, GFS_restart_var, Atm_block, GFS_control, Atmos%domain_for_read, & Atm(mygrid)%flagstruct%warm_start, ignore_rst_cksum) if(GFS_control%do_ca .and. Atm(mygrid)%flagstruct%warm_start)then @@ -1055,21 +1064,10 @@ subroutine atmos_model_end (Atmos) endif #endif - call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst) + call atmosphere_end (Atmos % Time, Atmos%grid, .false.) - if(restart_endfcst) then - call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, & - GFS_control, Atmos%domain) -! call write_stoch_restart_atm('RESTART/atm_stoch.res.nc') - endif if (GFS_Control%do_sppt .or. GFS_Control%do_shum .or. GFS_Control%do_skeb .or. & GFS_Control%lndp_type > 0 .or. GFS_Control%do_ca .or. GFS_Control%do_spp) then - if(restart_endfcst) then - call write_stoch_restart_atm('RESTART/atm_stoch.res.nc') - if (GFS_control%do_ca)then - call write_ca_restart() - endif - endif call stochastic_physics_wrapper_end(GFS_control) endif @@ -1099,9 +1097,15 @@ subroutine atmos_model_restart(Atmos, timestamp) type (atmos_data_type), intent(inout) :: Atmos character(len=*), intent(in) :: timestamp - call atmosphere_restart(timestamp) - call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, & - GFS_control, Atmos%domain, timestamp) + if (quilting_restart) then + call fv_sfc_restart_output(GFS_Data%Sfcprop, Atm_block, GFS_control) + call fv_phy_restart_output(GFS_restart_var, Atm_block) + call fv_dyn_restart_output(Atm(mygrid), timestamp) + else + call atmosphere_restart(timestamp) + call FV3GFS_restart_write (GFS_data, GFS_restart_var, Atm_block, & + GFS_control, Atmos%domain, timestamp) + endif if(GFS_control%do_ca)then call write_ca_restart(timestamp) endif diff --git a/ccpp/data/CCPP_typedefs.meta b/ccpp/data/CCPP_typedefs.meta index c7c5b4521..09c4feda5 100644 --- a/ccpp/data/CCPP_typedefs.meta +++ b/ccpp/data/CCPP_typedefs.meta @@ -739,20 +739,6 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys -[evbs] - standard_name = soil_upward_latent_heat_flux - long_name = soil upward latent heat flux - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys -[evcw] - standard_name = canopy_upward_latent_heat_flux - long_name = canopy upward latent heat flux - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys [pah] standard_name = total_precipitation_advected_heat long_name = precipitation advected heat - total @@ -1900,13 +1886,6 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys -[sbsno] - standard_name = snow_deposition_sublimation_upward_latent_heat_flux - long_name = latent heat flux from snow depo/subl - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys [scmpsw] standard_name = components_of_surface_downward_shortwave_fluxes long_name = derived type for special components of surface downward shortwave fluxes @@ -2107,13 +2086,6 @@ units = flag dimensions = () type = logical -[trans] - standard_name = transpiration_flux - long_name = total plant transpiration rate - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys [tseal] standard_name = surface_skin_temperature_for_nsst long_name = ocean surface skin temperature diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90 index f5059271b..56b01aa12 100644 --- a/ccpp/data/GFS_typedefs.F90 +++ b/ccpp/data/GFS_typedefs.F90 @@ -214,6 +214,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: lakefrac(:) => null() !< lake fraction [0:1] real (kind=kind_phys), pointer :: lakedepth(:) => null() !< lake depth [ m ] real (kind=kind_phys), pointer :: tsfc (:) => null() !< surface air temperature in K + real (kind=kind_phys), pointer :: vegtype_frac (:,:) => null() !< fractions [0:1] of veg. categories + real (kind=kind_phys), pointer :: soiltype_frac(:,:) => null() !< fractions [0:1] of soil categories !< [tsea in gbphys.f] real (kind=kind_phys), pointer :: tsfco (:) => null() !< sst in K real (kind=kind_phys), pointer :: tsfcl (:) => null() !< surface land temperature in K @@ -373,10 +375,13 @@ module GFS_typedefs real (kind=kind_phys), pointer :: clw_surf_ice(:) => null() !< RUC LSM: moist cloud water mixing ratio at surface over ice real (kind=kind_phys), pointer :: qwv_surf_land(:) => null() !< RUC LSM: water vapor mixing ratio at surface over land real (kind=kind_phys), pointer :: qwv_surf_ice(:) => null() !< RUC LSM: water vapor mixing ratio at surface over ice + real (kind=kind_phys), pointer :: rhofr(:) => null() !< RUC LSM: internal density of frozen precipitation real (kind=kind_phys), pointer :: tsnow_land(:) => null() !< RUC LSM: snow temperature at the bottom of the first snow layer over land real (kind=kind_phys), pointer :: tsnow_ice(:) => null() !< RUC LSM: snow temperature at the bottom of the first snow layer over ice real (kind=kind_phys), pointer :: snowfallac_land(:) => null() !< ruc lsm diagnostics over land real (kind=kind_phys), pointer :: snowfallac_ice(:) => null() !< ruc lsm diagnostics over ice + real (kind=kind_phys), pointer :: acsnow_land(:) => null() !< ruc lsm diagnostics over land + real (kind=kind_phys), pointer :: acsnow_ice(:) => null() !< ruc lsm diagnostics over ice ! MYNN surface layer real (kind=kind_phys), pointer :: ustm (:) => null() !u* including drag @@ -954,9 +959,11 @@ module GFS_typedefs !< ivegsrc = 3 => NLCD40 (40 category, NOAH WRFv4 only) !< ivegsrc = 4 => USGS-RUC (28 category, NOAH WRFv4 only) !< ivegsrc = 5 => MODI-RUC (21 category, NOAH WRFv4 only) + integer :: nvegcat !< nvegcat = 20 if ivegsrc = 1 integer :: isot !< isot = 0 => Zobler soil type ( 9 category) !< isot = 1 => STATSGO soil type (19 category, AKA 'STAS'(?)) !< isot = 2 => STAS-RUC soil type (19 category, NOAH WRFv4 only) + integer :: nsoilcat !< nsoilcat = 19 if isot = 1 integer :: kice !< number of layers in sice integer :: lsoil_lsm !< number of soil layers internal to land surface model integer :: lsnow_lsm !< maximum number of snow layers internal to land surface model @@ -991,6 +998,12 @@ module GFS_typedefs integer :: iopt_stc !snow/soil temperature time scheme (only layer 1) integer :: iopt_trs !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb inversed) + ! -- RUC LSM options + integer :: mosaic_lu=0 !< control for use of fractional landuse in RUC land surface model + integer :: mosaic_soil=0 !< control for use of fractional soil in RUC land surface model + integer :: isncond_opt=1 !< control for soil thermal conductivity option in RUC land surface model + integer :: isncovr_opt=1 !< control for snow cover fraction option in RUC land surface model + logical :: use_ufo !< flag for gcycle surface option ! GFDL Surface Layer options @@ -1206,6 +1219,10 @@ module GFS_typedefs !--- potential temperature definition in surface layer physics logical :: thsfc_loc !< flag for local vs. standard potential temperature +!--- flux method in 2-m diagnostics + logical :: diag_flux !< flag for flux method in 2-m diagnostics +!--- log method in 2-m diagnostics (for stable conditions) + logical :: diag_log !< flag for log method in 2-m diagnostics (for stable conditions) !--- vertical diffusion real(kind=kind_phys) :: xkzm_m !< [in] bkgd_vdif_m background vertical diffusion for momentum @@ -1264,6 +1281,10 @@ module GFS_typedefs ! 2 - scheme from Draper, JHM, 2021. real(kind=kind_phys) :: sppt_amp ! pjp cloud perturbations integer :: n_var_lndp + logical :: lndp_each_step ! flag to indicate that land perturbations are applied at every time step, + ! otherwise they are applied only + ! after gcycle is run + ! next two are duplicated here to support lndp_type=1. If delete that scheme, could remove from GFS defs? character(len=3) , pointer :: lndp_var_list(:) real(kind=kind_phys), pointer :: lndp_prt_list(:) @@ -1724,15 +1745,21 @@ module GFS_typedefs ! %upfx0 - clear sky upward lw flux at toa (w/m**2) ! Input/output - used by physics - real (kind=kind_phys), pointer :: srunoff(:) => null() !< surface water runoff (from lsm) - real (kind=kind_phys), pointer :: evbsa (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: evcwa (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: snohfa (:) => null() !< noah lsm diagnostics + real (kind=kind_phys), pointer :: srunoff(:) => null() !< accumulated surface storm runoff (from lsm) + real (kind=kind_phys), pointer :: evbsa (:) => null() !< accumulated direct evaporation + real (kind=kind_phys), pointer :: evcwa (:) => null() !< accumulated canopy evaporation + real (kind=kind_phys), pointer :: snohfa (:) => null() !< heat flux for phase change of snow (melting) + real (kind=kind_phys), pointer :: transa (:) => null() !< accumulated transpiration + real (kind=kind_phys), pointer :: sbsnoa (:) => null() !< accumulated snow sublimation + real (kind=kind_phys), pointer :: snowca (:) => null() !< snow cover + real (kind=kind_phys), pointer :: sbsno (:) => null() !< instantaneous snow sublimation + real (kind=kind_phys), pointer :: evbs(:) => null() !< instantaneous direct evaporation + real (kind=kind_phys), pointer :: trans (:) => null() !< instantaneous transpiration + real (kind=kind_phys), pointer :: evcw(:) => null() !< instantaneous canopy evaporation + real (kind=kind_phys), pointer :: snowmt_land(:) => null() !< ruc lsm diagnostics over land + real (kind=kind_phys), pointer :: snowmt_ice(:) => null() !< ruc lsm diagnostics over ice + real (kind=kind_phys), pointer :: soilm (:) => null() !< integrated soil moisture real (kind=kind_phys), pointer :: paha (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: transa (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: sbsnoa (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: snowca (:) => null() !< noah lsm diagnostics - real (kind=kind_phys), pointer :: soilm (:) => null() !< soil moisture real (kind=kind_phys), pointer :: tmpmin (:) => null() !< min temperature at 2m height (k) real (kind=kind_phys), pointer :: tmpmax (:) => null() !< max temperature at 2m height (k) real (kind=kind_phys), pointer :: dusfc (:) => null() !< u component of surface stress @@ -1783,7 +1810,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: tsnowp (:) => null() !< accumulated surface snowfall (m) real (kind=kind_phys), pointer :: tsnowpb(:) => null() !< accumulated surface snowfall in bucket (m) real (kind=kind_phys), pointer :: rhonewsn1(:) => null() !< precipitation ice density outside RUC LSM (kg/m3) - real (kind=kind_phys), pointer :: rhosnf(:) => null() !< precipitation ice density inside RUC LSM (kg/m3) !--- MYNN variables real (kind=kind_phys), pointer :: edmf_a (:,:) => null() ! @@ -2142,6 +2168,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model) allocate (Sfcprop%slmsk (IM)) allocate (Sfcprop%oceanfrac(IM)) allocate (Sfcprop%landfrac (IM)) + allocate (Sfcprop%vegtype_frac (IM,Model%nvegcat)) + allocate (Sfcprop%soiltype_frac(IM,Model%nsoilcat)) allocate (Sfcprop%lakefrac (IM)) allocate (Sfcprop%lakedepth(IM)) allocate (Sfcprop%tsfc (IM)) @@ -2171,10 +2199,14 @@ subroutine sfcprop_create (Sfcprop, IM, Model) allocate (Sfcprop%emis_lnd (IM)) allocate (Sfcprop%emis_ice (IM)) allocate (Sfcprop%emis_wat (IM)) + allocate (Sfcprop%acsnow_land (IM)) + allocate (Sfcprop%acsnow_ice (IM)) Sfcprop%slmsk = clear_val Sfcprop%oceanfrac = clear_val Sfcprop%landfrac = clear_val + Sfcprop%vegtype_frac = clear_val + Sfcprop%soiltype_frac = clear_val Sfcprop%lakefrac = clear_val Sfcprop%lakedepth = clear_val Sfcprop%tsfc = clear_val @@ -2204,6 +2236,9 @@ subroutine sfcprop_create (Sfcprop, IM, Model) Sfcprop%emis_lnd = clear_val Sfcprop%emis_ice = clear_val Sfcprop%emis_wat = clear_val + Sfcprop%acsnow_land = clear_val + Sfcprop%acsnow_ice = clear_val + !--- In (radiation only) allocate (Sfcprop%snoalb (IM)) @@ -2487,10 +2522,13 @@ subroutine sfcprop_create (Sfcprop, IM, Model) allocate (Sfcprop%clw_surf_ice (IM)) allocate (Sfcprop%qwv_surf_land (IM)) allocate (Sfcprop%qwv_surf_ice (IM)) + allocate (Sfcprop%rhofr (IM)) allocate (Sfcprop%tsnow_land (IM)) allocate (Sfcprop%tsnow_ice (IM)) allocate (Sfcprop%snowfallac_land (IM)) allocate (Sfcprop%snowfallac_ice (IM)) + allocate (Sfcprop%acsnow_land (IM)) + allocate (Sfcprop%acsnow_ice (IM)) ! Sfcprop%wetness = clear_val Sfcprop%sh2o = clear_val @@ -2502,10 +2540,13 @@ subroutine sfcprop_create (Sfcprop, IM, Model) Sfcprop%qwv_surf_land = clear_val Sfcprop%qwv_surf_ice = clear_val Sfcprop%flag_frsoil = clear_val + Sfcprop%rhofr = -1.e3 Sfcprop%tsnow_land = clear_val Sfcprop%tsnow_ice = clear_val Sfcprop%snowfallac_land = clear_val Sfcprop%snowfallac_ice = clear_val + Sfcprop%acsnow_land = clear_val + Sfcprop%acsnow_ice = clear_val ! if (Model%rdlai) then allocate (Sfcprop%xlaixy (IM)) @@ -3241,8 +3282,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: ivegsrc = 2 !< ivegsrc = 0 => USGS, !< ivegsrc = 1 => IGBP (20 category) !< ivegsrc = 2 => UMD (13 category) + integer :: nvegcat = 20 !< number of veg. categories depending on ivegsrc integer :: isot = 0 !< isot = 0 => Zobler soil type ( 9 category) !< isot = 1 => STATSGO soil type (19 category) + integer :: nsoilcat = 16 !< number of soil categories depending on isot ! -- to use Noah MP, lsm needs to be set to 2 and both ivegsrc and isot are set ! to 1 - MODIS IGBP and STATSGO - the defaults are the same as in the ! scripts;change from namelist @@ -3261,6 +3304,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iopt_stc = 1 !snow/soil temperature time scheme (only layer 1) integer :: iopt_trs = 2 !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb reversed) + integer :: mosaic_lu = 0 ! 1 - used of fractional landuse in RUC lsm + integer :: mosaic_soil = 0 ! 1 - used of fractional soil in RUC lsm + integer :: isncond_opt = 1 ! 2 - Sturm (1997) + integer :: isncovr_opt = 1 ! 2 - Niu-Yang (2007), 3-updated Niu-Yang similar to Noah MP + logical :: use_ufo = .false. !< flag for gcycle surface option logical :: lcurr_sf = .false. !< flag for taking ocean currents into account in GFDL surface layer @@ -3469,6 +3517,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- potential temperature definition in surface layer physics logical :: thsfc_loc = .true. !< flag for local vs. standard potential temperature +!--- flux method in 2-m diagnostics + logical :: diag_flux = .false. !< flag for flux method in 2-m diagnostics +!--- flux method in 2-m diagnostics (for stable conditions) + logical :: diag_log = .false. !< flag for log method in 2-m diagnostics (for stable conditions) !<.true. means use local (gridpoint) surface pressure to define potential temperature !< this is the current GFS physics approach !<.false. means use reference pressure of 1000 hPa to define potential temperature @@ -3549,6 +3601,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: skeb_npass = 11 integer :: lndp_type = 0 integer :: n_var_lndp = 0 + logical :: lndp_each_step = .false. integer :: n_var_spp = 0 integer :: spp_pbl = 0 integer :: spp_sfc = 0 @@ -3641,11 +3694,13 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- land/surface model control lsm, lsoil, lsoil_lsm, lsnow_lsm, kice, rdlai, & nmtvr, ivegsrc, use_ufo, iopt_thcnd, ua_phys, usemonalb, & - aoasis, fasdas,exticeden, & + aoasis, fasdas, exticeden, nvegcat, nsoilcat, & ! Noah MP options iopt_dveg,iopt_crs,iopt_btr,iopt_run,iopt_sfc, iopt_frz, & iopt_inf, iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc, & iopt_trs, & + ! RUC lsm options + mosaic_lu, mosaic_soil, isncond_opt, isncovr_opt, & ! GFDL surface layer options lcurr_sf, pert_cd, ntsflg, sfenth, & !--- lake model control @@ -3677,7 +3732,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & dlqf, rbcr, shoc_parm, psauras, prauras, wminras, & do_sppt, do_shum, do_skeb, & do_spp, n_var_spp, & - lndp_type, n_var_lndp, & + lndp_type, n_var_lndp, lndp_each_step, & pert_mp,pert_clds,pert_radtend, & !--- Rayleigh friction prslrd0, ral_ts, ldiag_ugwp, do_ugwp, do_tofd, & @@ -3697,6 +3752,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & sfc_z0_type, & !--- switch beteeen local and standard potential temperature thsfc_loc, & + !--- switches in 2-m diagnostics + 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, & @@ -4340,7 +4397,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%aoasis = aoasis Model%fasdas = fasdas Model%ivegsrc = ivegsrc + Model%nvegcat = nvegcat Model%isot = isot + Model%nsoilcat = nsoilcat Model%use_ufo = use_ufo Model%exticeden = exticeden if (Model%exticeden .and. & @@ -4380,6 +4439,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%iopt_stc = iopt_stc Model%iopt_trs = iopt_trs +! RUC lsm options + Model%mosaic_lu = mosaic_lu + Model%mosaic_soil = mosaic_soil + Model%isncond_opt = isncond_opt + Model%isncovr_opt = isncovr_opt + !--- tuning parameters for physical parameterizations Model%ras = ras Model%flipv = flipv @@ -4568,6 +4633,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- potential temperature reference in sfc layer Model%thsfc_loc = thsfc_loc +!--- flux method in 2-m diagnostics + Model%diag_flux = diag_flux +!--- flux method in 2-m diagnostics (for stable conditions) + Model%diag_log = diag_log !--- vertical diffusion Model%xkzm_m = xkzm_m @@ -4601,6 +4670,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- stochastic surface perturbation options Model%lndp_type = lndp_type Model%n_var_lndp = n_var_lndp + Model%lndp_each_step = lndp_each_step Model%do_spp = do_spp Model%n_var_spp = n_var_spp @@ -5295,6 +5365,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & print *,'iopt_trs = ', Model%iopt_trs elseif (Model%lsm == Model%lsm_ruc) then print *,' RUC Land Surface Model used' + print *, 'The Physics options are' + print *,' mosaic_lu = ',mosaic_lu + print *,' mosaic_soil = ',mosaic_soil + print *,' isncond_opt = ',isncond_opt + print *,' isncovr_opt = ',isncovr_opt else print *,' Unsupported LSM type - job aborted - lsm=',Model%lsm stop @@ -6160,7 +6235,9 @@ subroutine control_print(Model) print *, ' shape(pores) : ', shape(Model%pores) print *, ' shape(resid) : ', shape(Model%resid) print *, ' ivegsrc : ', Model%ivegsrc + print *, ' nvegcat : ', Model%nvegcat print *, ' isot : ', Model%isot + print *, ' nsoilcat : ', Model%nsoilcat if (Model%lsm == Model%lsm_noahmp) then print *, ' Noah MP LSM is used, the options are' @@ -6177,6 +6254,13 @@ subroutine control_print(Model) print *, ' iopt_tbot : ', Model%iopt_tbot print *, ' iopt_stc : ', Model%iopt_stc print *, ' iopt_trs : ', Model%iopt_trs + elseif (Model%lsm == Model%lsm_ruc) then + print *,' RUC Land Surface Model used' + print *, 'The Physics options are' + print *,' mosaic_lu = ',Model%mosaic_lu + print *,' mosaic_soil = ',Model%mosaic_soil + print *,' isncond_opt = ',Model%isncond_opt + print *,' isncovr_opt = ',Model%isncovr_opt endif print *, ' use_ufo : ', Model%use_ufo print *, ' lcurr_sf : ', Model%lcurr_sf @@ -6246,6 +6330,8 @@ subroutine control_print(Model) print *, ' rbcr : ', Model%rbcr print *, ' do_mynnedmf : ', Model%do_mynnedmf print *, ' do_mynnsfclay : ', Model%do_mynnsfclay + print *, ' diag_flux : ', Model%diag_flux + print *, ' diag_log : ', Model%diag_log print *, ' do_myjsfc : ', Model%do_myjsfc print *, ' do_myjpbl : ', Model%do_myjpbl print *, ' do_ugwp : ', Model%do_ugwp @@ -6326,6 +6412,7 @@ subroutine control_print(Model) print *, ' do_skeb : ', Model%do_skeb print *, ' lndp_type : ', Model%lndp_type print *, ' n_var_lndp : ', Model%n_var_lndp + print *, ' lndp_each_step : ', Model%lndp_each_step print *, ' do_spp : ', Model%do_spp print *, ' n_var_spp : ', Model%n_var_spp print *, ' ' @@ -7053,6 +7140,12 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%transa (IM)) allocate (Diag%sbsnoa (IM)) allocate (Diag%snowca (IM)) + allocate (Diag%evbs (IM)) + allocate (Diag%evcw (IM)) + allocate (Diag%sbsno (IM)) + allocate (Diag%trans (IM)) + allocate (Diag%snowmt_land (IM)) + allocate (Diag%snowmt_ice (IM)) allocate (Diag%soilm (IM)) allocate (Diag%tmpmin (IM)) allocate (Diag%tmpmax (IM)) @@ -7122,7 +7215,6 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%smcwlt2 (IM)) allocate (Diag%smcref2 (IM)) allocate (Diag%rhonewsn1 (IM)) - allocate (Diag%rhosnf (IM)) allocate (Diag%frzr (IM)) allocate (Diag%frzrb (IM)) allocate (Diag%frozr (IM)) @@ -7370,8 +7462,14 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%evcwa = zero Diag%snohfa = zero Diag%transa = zero - Diag%sbsnoa = zero Diag%snowca = zero + Diag%sbsnoa = zero + Diag%sbsno = zero + Diag%evbs = zero + Diag%evcw = zero + Diag%trans = zero + Diag%snowmt_land= zero + Diag%snowmt_ice = zero Diag%soilm = zero Diag%tmpmin = Model%huge Diag%tmpmax = zero @@ -7589,7 +7687,6 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%rh02min = 999. Diag%pratemax = 0. Diag%rhonewsn1 = 200. - Diag%rhosnf = -1.e3 set_totprcp = .false. if (present(linit) ) set_totprcp = linit if (present(iauwindow_center) ) set_totprcp = iauwindow_center diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta index 9dfbf378f..cdca4cacf 100644 --- a/ccpp/data/GFS_typedefs.meta +++ b/ccpp/data/GFS_typedefs.meta @@ -613,6 +613,20 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys +[vegtype_frac] + standard_name = fraction_of_vegetation_category + long_name = fraction of horizontal grid area occupied by given vegetation category + units = frac + dimensions = (horizontal_loop_extent,number_of_vegetation_categories) + type = real + kind = kind_phys +[soiltype_frac] + standard_name = fraction_of_soil_category + long_name = fraction of horizontal grid area occupied by given soil category + units = frac + dimensions = (horizontal_loop_extent,number_of_soil_categories) + type = real + kind = kind_phys [lakefrac] standard_name = lake_area_fraction long_name = fraction of horizontal grid area occupied by lake @@ -1635,6 +1649,14 @@ type = real kind = kind_phys active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme) +[rhofr] + standard_name = lsm_internal_surface_frozen_precipitation_density + long_name = density of frozen precipitation + units = kg m-3 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme) [tsnow_land] standard_name = temperature_in_surface_snow_at_surface_adjacent_layer_over_land long_name = snow temperature at the bottom of the first snow layer over land @@ -1652,21 +1674,35 @@ kind = kind_phys active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme) [snowfallac_land] - standard_name = surface_snow_amount_over_land - long_name = run-total snow accumulation on the ground + standard_name = surface_snow_amount_assuming_variable_snow_density_over_land + long_name = run-total snow accumulation on the ground with variable snow density over land units = kg m-2 dimensions = (horizontal_loop_extent) type = real kind = kind_phys active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme) +[acsnow_land] + standard_name = surface_snow_lwe_thickness_amount_over_land + long_name = run-total snowfall water equivalent over land + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys [snowfallac_ice] - standard_name = surface_snow_amount_over_ice - long_name = run-total snow accumulation on the ice + standard_name = surface_snow_amount_assuming_variable_snow_density_over_ice + long_name = run-total snow accumulation on the ground with variable snow density over ice units = kg m-2 dimensions = (horizontal_loop_extent) type = real kind = kind_phys active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme) +[acsnow_ice] + standard_name = surface_snow_lwe_thickness_amount_over_ice + long_name = run-total snowfall water equivalent over ice + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys [ustm] standard_name = surface_friction_velocity_for_momentum long_name = friction velocity isolated for momentum only @@ -4266,6 +4302,30 @@ units = flag dimensions = () type = integer +[mosaic_lu] + standard_name = control_for_fractional_landuse_in_ruc_land_surface_scheme + long_name = control for use of fractional landuse info in RUC land surface model + units = flag + dimensions = () + type = integer +[mosaic_soil] + standard_name = control_for_fractional_soil_in_ruc_land_surface_scheme + long_name = control for use of fractional soil info in RUC land surface model + units = flag + dimensions = () + type = integer +[isncond_opt] + standard_name = control_for_soil_thermal_conductivity_option_in_ruc_lsm + long_name = control for soil thermal conductivity option in RUC land surface model + units = flag + dimensions = () + type = integer +[isncovr_opt] + standard_name = control_for_snow_cover_fraction_option_in_ruc_lsm + long_name = control for snow cover fraction option in RUC land surface model + units = flag + dimensions = () + type = integer [kice] standard_name = vertical_dimension_of_sea_ice long_name = vertical loop extent for ice levels, start at 1 @@ -4373,12 +4433,24 @@ units = index dimensions = () type = integer +[nvegcat] + standard_name = number_of_vegetation_categories + long_name = number of vegetation categories + units = count + dimensions = () + type = integer [isot] standard_name = control_for_soil_type_dataset long_name = soil type dataset choice units = index dimensions = () type = integer +[nsoilcat] + standard_name = number_of_soil_categories + long_name = number of soil categories + units = count + dimensions = () + type = integer [iopt_thcnd] standard_name = control_for_land_surface_scheme_thermal_conductivity_option long_name = choice for thermal conductivity option (see module_sf_noahlsm) @@ -4646,6 +4718,18 @@ units = flag dimensions = () type = logical +[diag_flux] + standard_name = flag_for_flux_method_in_2m_diagnostics + long_name = flag for flux method in 2-m diagnostics + units = flag + dimensions = () + type = logical +[diag_log] + standard_name = flag_for_log_method_in_2m_diagnostics + long_name = flag for log method in 2-m diagnostics + units = flag + dimensions = () + type = logical [hybedmf] standard_name = flag_for_hybrid_edmf_pbl_scheme long_name = flag for hybrid edmf pbl scheme (moninedmf) @@ -7595,6 +7679,34 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys +[sbsno] + standard_name = snow_deposition_sublimation_upward_latent_heat_flux + long_name = latent heat flux from snow depo/subl + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys +[evbs] + standard_name = soil_upward_latent_heat_flux + long_name = soil upward latent heat flux + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys +[evcw] + standard_name = canopy_upward_latent_heat_flux + long_name = canopy upward latent heat flux + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys +[trans] + standard_name = transpiration_flux + long_name = total plant transpiration rate + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys [soilm] standard_name = soil_moisture_content long_name = soil moisture @@ -7602,6 +7714,20 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys +[snowmt_land] + standard_name = surface_snow_melt_over_land + long_name = snow melt during timestep over land + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys +[snowmt_ice] + standard_name = surface_snow_melt_over_ice + long_name = snow melt during timestep over ice + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys [tmpmin] standard_name = minimum_temperature_at_2m long_name = min temperature at 2m height @@ -7749,13 +7875,6 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys -[rhosnf] - standard_name = lsm_internal_surface_frozen_precipitation_density - long_name = density of frozen precipitation - units = kg m-3 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys [train] standard_name = accumulated_change_of_air_temperature_due_to_FA_scheme long_name = accumulated change of air temperature due to FA MP scheme diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90 index 5d693296e..83d20ea59 100644 --- a/ccpp/driver/GFS_diagnostics.F90 +++ b/ccpp/driver/GFS_diagnostics.F90 @@ -1060,7 +1060,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'ssrun_acc' - ExtDiag(idx)%desc = 'surface storm water runoff - GFS lsm' + ExtDiag(idx)%desc = 'Accumulated surface storm water runoff' ExtDiag(idx)%unit = 'kg/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' allocate (ExtDiag(idx)%data(nblks)) @@ -1071,7 +1071,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'evbs_ave' - ExtDiag(idx)%desc = 'Direct Evaporation from Bare Soil - GFS lsm' + ExtDiag(idx)%desc = 'Direct Evaporation from Bare Soil' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1083,7 +1083,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'evcw_ave' - ExtDiag(idx)%desc = 'Canopy water evaporation - GFS lsm' + ExtDiag(idx)%desc = 'Canopy water evaporation' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1095,7 +1095,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'snohf' - ExtDiag(idx)%desc = 'Snow Phase Change Heat Flux - GFS lsm' + ExtDiag(idx)%desc = 'Snow Phase Change Heat Flux' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1108,7 +1108,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'pah_ave' - ExtDiag(idx)%desc = ' Total Precipitation Advected Heat - GFS lsm' + ExtDiag(idx)%desc = ' Total Precipitation Advected Heat' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1121,7 +1121,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'trans_ave' - ExtDiag(idx)%desc = 'transpiration - GFS lsm' + ExtDiag(idx)%desc = 'transpiration' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1133,7 +1133,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'sbsno_ave' - ExtDiag(idx)%desc = 'Sublimation (evaporation from snow) - GFS lsm' + ExtDiag(idx)%desc = 'Sublimation (evaporation from snow)' ExtDiag(idx)%unit = 'W/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%time_avg = .TRUE. @@ -1155,6 +1155,17 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowca(:) enddo + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'snowc' + ExtDiag(idx)%desc = 'snow cover ' + ExtDiag(idx)%unit = 'fraction' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%sncovr(:) + enddo + idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'soilm' @@ -1295,7 +1306,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'dlwsfc' - ExtDiag(idx)%desc = 'time accumulated downward lw flux at surface- GFS physics' + ExtDiag(idx)%desc = 'time accumulated downward lw flux at surface' ExtDiag(idx)%unit = 'w/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%intpl_method = 'bilinear' @@ -1307,7 +1318,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'ulwsfc' - ExtDiag(idx)%desc = 'time accumulated upward lw flux at surface- GFS physics' + ExtDiag(idx)%desc = 'time accumulated upward lw flux at surface' ExtDiag(idx)%unit = 'w/m**2' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%intpl_method = 'bilinear' @@ -2192,7 +2203,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'cpofp' ExtDiag(idx)%desc = 'Percent frozen precipitation' - ExtDiag(idx)%unit = '%' + ExtDiag(idx)%unit = 'fraction' ExtDiag(idx)%mod_name = 'gfs_phys' ExtDiag(idx)%intpl_method = 'bilinear' allocate (ExtDiag(idx)%data(nblks)) @@ -3290,7 +3301,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'canopy' ExtDiag(idx)%desc = 'canopy water (cnwat in gfs data)' - ExtDiag(idx)%unit = 'XXX' + ExtDiag(idx)%unit = 'mm' ExtDiag(idx)%mod_name = 'gfs_sfc' allocate (ExtDiag(idx)%data(nblks)) do nb = 1,nblks @@ -3389,7 +3400,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'hice' ExtDiag(idx)%desc = 'sea ice thickness (icetk in gfs_data)' - ExtDiag(idx)%unit = 'XXX' + ExtDiag(idx)%unit = 'm' ExtDiag(idx)%mod_name = 'gfs_sfc' allocate (ExtDiag(idx)%data(nblks)) do nb = 1,nblks @@ -3441,7 +3452,62 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowd(:) enddo + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'sbsno' + ExtDiag(idx)%desc = 'instantaneous sublimation (evaporation from snow)' + ExtDiag(idx)%unit = 'W/m**2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sbsno(:) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'evbs' + ExtDiag(idx)%desc = 'instantaneous direct evaporation over land' + ExtDiag(idx)%unit = 'W m-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%evbs(:) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'evcw' + ExtDiag(idx)%desc = 'instantaneous canopy evaporation' + ExtDiag(idx)%unit = 'W m-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%evcw(:) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'trans' + ExtDiag(idx)%desc = 'instantaneous transpiration' + ExtDiag(idx)%unit = 'W m-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%trans(:) + enddo + if (Model%lsm == Model%lsm_ruc) then + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'rhofr' + ExtDiag(idx)%desc = 'density of frozen precipitation' + ExtDiag(idx)%unit = 'kg m-3' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%rhofr(:) + enddo + idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'snowfall_acc_land' @@ -3453,6 +3519,28 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowfallac_land(:) enddo + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'acsnow_land' + ExtDiag(idx)%desc = 'total accumulated SWE of frozen precipitation over land' + ExtDiag(idx)%unit = 'kg m-2' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%acsnow_land(:) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'snowmt_land' + ExtDiag(idx)%desc = 'accumulated snow melt over land' + ExtDiag(idx)%unit = 'kg m-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowmt_land(:) + enddo + idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'snowfall_acc_ice' @@ -3463,7 +3551,29 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowfallac_ice(:) enddo - endif + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'acsnow_ice' + ExtDiag(idx)%desc = 'total accumulated SWE of frozen precipitation over ice' + ExtDiag(idx)%unit = 'kg m-2' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%acsnow_ice(:) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'snowmt_ice' + ExtDiag(idx)%desc = 'accumulated snow melt over ice' + ExtDiag(idx)%unit = 'kg m-2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowmt_ice(:) + enddo + endif ! RUC lsm idx = idx + 1 ExtDiag(idx)%axes = 2 @@ -3575,8 +3685,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'tprcp' - ExtDiag(idx)%desc = 'total precipitation' - ExtDiag(idx)%unit = 'kg/m**2' + ExtDiag(idx)%desc = 'total time-step precipitation' + ExtDiag(idx)%unit = 'm' ExtDiag(idx)%mod_name = 'gfs_sfc' allocate (ExtDiag(idx)%data(nblks)) do nb = 1,nblks @@ -3586,7 +3696,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'vtype' - ExtDiag(idx)%desc = 'vegetation type in integer 1-13' + ExtDiag(idx)%desc = 'vegetation type in integer' ExtDiag(idx)%unit = 'number' ExtDiag(idx)%mod_name = 'gfs_sfc' allocate (ExtDiag(idx)%data(nblks)) @@ -3652,6 +3762,18 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%vfrac(:) enddo + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'wetness' + ExtDiag(idx)%desc = 'soil moisture availability in top soil layer' + ExtDiag(idx)%unit = 'fraction' + ExtDiag(idx)%mod_name = 'gfs_sfc' + ExtDiag(idx)%cnvfac = cn_100 + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%wetness(:) + enddo + if (Model%rdlai) then idx = idx + 1 ExtDiag(idx)%axes = 2 @@ -3665,6 +3787,34 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop enddo endif + do num = 1,Model%nvegcat + write (xtra,'(i2)') num + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'vfrac_'//trim(xtra) + ExtDiag(idx)%desc = 'fraction of vegetation category'//trim(xtra) + ExtDiag(idx)%unit = 'frac' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%vegtype_frac(:,num) + enddo + enddo + + do num = 1,Model%nsoilcat + write (xtra,'(i2)') num + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = 'sfrac_'//trim(xtra) + ExtDiag(idx)%desc = 'fraction of soil category'//trim(xtra) + ExtDiag(idx)%unit = 'frac' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%soiltype_frac(:,num) + enddo + enddo + if (Model%lsm == Model%lsm_ruc) then do num = 1,Model%lsoil_lsm write (xtra,'(i1)') num diff --git a/ccpp/physics b/ccpp/physics index ec115b818..b7dc3ec0c 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit ec115b81804073e3c4054de0ed7bfa6a1bca678a +Subproject commit b7dc3ec0c966ff76d8f40cce0eb4977bebea283b diff --git a/fv3_cap.F90 b/fv3_cap.F90 index 731d0e239..cdf53561b 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -27,7 +27,7 @@ module fv3gfs_cap_mod label_Finalize, & NUOPC_ModelGet ! - use module_fv3_config, only: quilting, output_fh, & + use module_fv3_config, only: quilting, quilting_restart, output_fh, & nfhout, nfhout_hf, nsout, dt_atmos, & calendar, cpl_grid_id, & cplprint_flag,output_1st_tstep_rst, & @@ -58,6 +58,7 @@ module fv3gfs_cap_mod type(ESMF_GridComp) :: fcstComp type(ESMF_State) :: fcstState type(ESMF_FieldBundle), allocatable :: fcstFB(:) + integer,dimension(:), allocatable :: fcstPetList integer, save :: FBCount type(ESMF_GridComp), allocatable :: wrtComp(:) @@ -73,6 +74,7 @@ module fv3gfs_cap_mod integer :: mype = -1 integer :: dbug = 0 + integer :: frestart(999) = -1 !----------------------------------------------------------------------- @@ -187,8 +189,8 @@ subroutine InitializeAdvertise(gcomp, rc) real :: nfhmax real :: output_startfh, outputfh, outputfh2(2) logical :: loutput_fh, lfreq - character(ESMF_MAXSTR) :: name - integer,dimension(:), allocatable :: petList, fcstPetList, originPetList, targetPetList + character(ESMF_MAXSTR) :: gc_name, fb_name + integer,dimension(:), allocatable :: petList, originPetList, targetPetList character(len=esmf_maxstr),allocatable :: fcstItemNameList(:) type(ESMF_StateItem_Flag), allocatable :: fcstItemTypeList(:) character(20) :: cwrtcomp @@ -205,9 +207,10 @@ subroutine InitializeAdvertise(gcomp, rc) character(len=*),parameter :: subname='(fv3_cap:InitializeAdvertise)' real(kind=8) :: MPI_Wtime, timeis, timerhs - integer :: wrttasks_per_group_from_parent, wrtLocalPet + integer :: wrttasks_per_group_from_parent, wrtLocalPet, num_threads character(len=64) :: rh_filename logical :: use_saved_routehandles, rh_file_exist + logical :: fieldbundle_is_restart = .false. ! !------------------------------------------------------------------------ @@ -215,12 +218,18 @@ subroutine InitializeAdvertise(gcomp, rc) rc = ESMF_SUCCESS timeis = MPI_Wtime() - call ESMF_GridCompGet(gcomp,name=name,vm=vm,rc=rc) + call ESMF_GridCompGet(gcomp, name=gc_name, vm=vm,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_VMGet(vm, petCount=petcount, localpet=mype, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! num_threads is needed to compute actual wrttasks_per_group_from_parent + call ESMF_InfoGetFromHost(gcomp, info=info, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_InfoGet(info, key="/NUOPC/Hint/PePerPet/MaxCount", value=num_threads, default=1, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! query for importState and exportState call NUOPC_ModelGet(gcomp, driverClock=clock, importState=importState, exportState=exportState, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -269,6 +278,12 @@ subroutine InitializeAdvertise(gcomp, rc) label ='quilting:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_ConfigGetAttribute(config=CF,value=quilting_restart, & + default=.false., label ='quilting_restart:',rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (.not.quilting) quilting_restart = .false. + call ESMF_ConfigGetAttribute(config=CF,value=output_1st_tstep_rst, & default=.false., label ='output_1st_tstep_rst:',rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -351,6 +366,7 @@ subroutine InitializeAdvertise(gcomp, rc) ! create fcst grid component if( quilting ) then + wrttasks_per_group_from_parent = wrttasks_per_group_from_parent * num_threads num_pes_fcst = petcount - write_groups * wrttasks_per_group_from_parent else num_pes_fcst = petcount @@ -396,7 +412,7 @@ subroutine InitializeAdvertise(gcomp, rc) ! determine number elements in fcstState call ESMF_StateGet(fcstState, itemCount=FBCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - if(mype == 0) print *,'af fcstCom FBCount= ',FBcount + if(mype == 0) print *,'fv3_cap: field bundles in fcstComp export state, FBCount= ',FBcount ! ! set start time for output output_startfh = 0. @@ -461,6 +477,9 @@ subroutine InitializeAdvertise(gcomp, rc) endif call ESMF_AttributeGet(fcstFB(i), convention="NetCDF", purpose="FV3", name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_AttributeGet(fcstFB(i), convention="NetCDF", purpose="FV3-nooutput", name="frestart", valueList=frestart, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + is_moving_fb(i) = is_moving(grid_id) enddo ! @@ -540,8 +559,7 @@ subroutine InitializeAdvertise(gcomp, rc) if(mype==0) print *,'af wrtState reconcile, FBcount=',FBcount - call ESMF_AttributeCopy(fcstState, wrtState(i), & - attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) + call ESMF_AttributeCopy(fcstState, wrtState(i), attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! deal with GridTransfer if needed @@ -596,8 +614,8 @@ subroutine InitializeAdvertise(gcomp, rc) ! loop over all FieldBundle in the states, for moving nests handle GridTransfer do j=1, FBcount if (is_moving_fb(j)) then - ! access the fcst (provider) Grid - call ESMF_FieldBundleGet(fcstFB(j), grid=providerGrid, rc=rc) + ! access the fcst (provider) Grid and fieldbundle name + call ESMF_FieldBundleGet(fcstFB(j), grid=providerGrid, name=fb_name, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! access the mirror FieldBundle on the wrtComp call ESMF_StateGet(wrtState(i), itemName="mirror_"//trim(fcstItemNameList(j)), fieldbundle=mirrorFB, rc=rc) @@ -683,22 +701,30 @@ subroutine InitializeAdvertise(gcomp, rc) if(mype == 0) print *,'af get wrtfb=',"output_"//trim(fcstItemNameList(j)),' rc=',rc if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - ! determine regridmethod - if (index(fcstItemNameList(j),"_bilinear") >0 ) then - regridmethod = ESMF_REGRIDMETHOD_BILINEAR - else if (index(fcstItemNameList(j),"_patch") >0) then - regridmethod = ESMF_REGRIDMETHOD_PATCH - else if (index(fcstItemNameList(j),"_nearest_stod") >0) then - regridmethod = ESMF_REGRIDMETHOD_NEAREST_STOD - else if (index(fcstItemNameList(j),"_nearest_dtos") >0) then - regridmethod = ESMF_REGRIDMETHOD_NEAREST_DTOS - else if (index(fcstItemNameList(j),"_conserve") >0) then - regridmethod = ESMF_REGRIDMETHOD_CONSERVE + fieldbundle_is_restart = .false. + if (fcstItemNameList(j)(1:8) == "restart_") then + ! restart output forecast bundles, no need to set regridmethod + ! Redist will be used instead of Regrid + fieldbundle_is_restart = .true. else - call ESMF_LogSetError(ESMF_RC_ARG_BAD, & - msg="Unable to determine regrid method.", & - line=__LINE__, file=__FILE__, rcToReturn=rc) - return + ! history output forecast bundles + ! determine regridmethod + if (index(fcstItemNameList(j),"_bilinear") >0 ) then + regridmethod = ESMF_REGRIDMETHOD_BILINEAR + else if (index(fcstItemNameList(j),"_patch") >0) then + regridmethod = ESMF_REGRIDMETHOD_PATCH + else if (index(fcstItemNameList(j),"_nearest_stod") >0) then + regridmethod = ESMF_REGRIDMETHOD_NEAREST_STOD + else if (index(fcstItemNameList(j),"_nearest_dtos") >0) then + regridmethod = ESMF_REGRIDMETHOD_NEAREST_DTOS + else if (index(fcstItemNameList(j),"_conserve") >0) then + regridmethod = ESMF_REGRIDMETHOD_CONSERVE + else + call ESMF_LogSetError(ESMF_RC_ARG_BAD, & + msg="Unable to determine regrid method.", & + line=__LINE__, file=__FILE__, rcToReturn=rc) + return + endif endif call ESMF_LogWrite('bf FieldBundleRegridStore', ESMF_LOGMSG_INFO, rc=rc) @@ -716,21 +742,33 @@ subroutine InitializeAdvertise(gcomp, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return else ! this is a Store() for the first wrtComp -> must do the Store() - call ESMF_TraceRegionEnter("ESMF_FieldBundleRegridStore()", rc=rc) - call ESMF_FieldBundleRegridStore(fcstFB(j), wrtFB(j,1), & - regridMethod=regridmethod, routehandle=routehandle(j,1), & - unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & - srcTermProcessing=isrcTermProcessing, rc=rc) - - ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - if (rc /= ESMF_SUCCESS) then - write(0,*)'fv3_cap.F90:InitializeAdvertise error in ESMF_FieldBundleRegridStore' - call ESMF_LogWrite('fv3_cap.F90:InitializeAdvertise error in ESMF_FieldBundleRegridStore', ESMF_LOGMSG_ERROR, rc=rc) - call ESMF_Finalize(endflag=ESMF_END_ABORT) + if (fieldbundle_is_restart) then + call ESMF_TraceRegionEnter("ESMF_FieldBundleRedistStore()", rc=rc) + call ESMF_FieldBundleRedistStore(fcstFB(j), wrtFB(j,1), & + routehandle=routehandle(j,1), & + rc=rc) + if (rc /= ESMF_SUCCESS) then + call ESMF_LogWrite('fv3_cap.F90:InitializeAdvertise error in ESMF_FieldBundleRedistStore', ESMF_LOGMSG_ERROR, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + call ESMF_TraceRegionExit("ESMF_FieldBundleRedistStore()", rc=rc) + call ESMF_LogWrite('af FieldBundleRedistStore', ESMF_LOGMSG_INFO, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + call ESMF_TraceRegionEnter("ESMF_FieldBundleRegridStore()", rc=rc) + call ESMF_FieldBundleRegridStore(fcstFB(j), wrtFB(j,1), & + regridMethod=regridmethod, routehandle=routehandle(j,1), & + unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & + srcTermProcessing=isrcTermProcessing, rc=rc) + if (rc /= ESMF_SUCCESS) then + call ESMF_LogWrite('fv3_cap.F90:InitializeAdvertise error in ESMF_FieldBundleRegridStore', ESMF_LOGMSG_ERROR, rc=rc) + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + call ESMF_TraceRegionExit("ESMF_FieldBundleRegridStore()", rc=rc) + call ESMF_LogWrite('af FieldBundleRegridStore', ESMF_LOGMSG_INFO, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif - call ESMF_TraceRegionExit("ESMF_FieldBundleRegridStore()", rc=rc) - call ESMF_LogWrite('af FieldBundleRegridStore', ESMF_LOGMSG_INFO, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (use_saved_routehandles) then call ESMF_RouteHandleWrite(routehandle(j,1), fileName=trim(rh_filename), rc=rc) @@ -916,7 +954,16 @@ subroutine InitializeAdvertise(gcomp, rc) endif ! end loutput_fh endif if(mype==0) print *,'output_fh=',output_fh(1:size(output_fh)),'lflname_fulltime=',lflname_fulltime -! + + if ( quilting ) then + do i=1, write_groups + call ESMF_InfoGetFromHost(wrtState(i), info=info, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_InfoSet(info, key="output_fh", values=output_fh, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + enddo + endif + ! --- advertise Fields in importState and exportState ------------------- ! call fcst Initialize (advertise phase) @@ -1085,8 +1132,8 @@ subroutine ModelAdvance_phase2(gcomp, rc) na = nint(time_elapsed/timeStep) call ESMF_TimeIntervalGet(time_elapsed, s=nfseconds, rc=rc) - output: if (ANY(nint(output_fh(:)*3600.0) == nfseconds)) then -! + output: if (ANY(nint(output_fh(:)*3600.0) == nfseconds) .or. ANY(frestart(:) == nfseconds)) then + if (mype == 0 .or. mype == lead_wrttask(1)) print *,' aft fcst run output time=',nfseconds, & 'FBcount=',FBcount,'na=',na diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index eb4de13e0..a84eabcfb 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -72,7 +72,7 @@ module FV3GFS_io_mod character(len=32), allocatable, dimension(:) :: oro_name2, sfc_name2, sfc_name3 real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: oro_var2, sfc_var2, phy_var2, sfc_var3ice character(len=32), allocatable, dimension(:) :: oro_ls_ss_name - real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: oro_ls_var, oro_ss_var + real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: oro_ls_var, oro_ss_var, oro_var3v, oro_var3s real(kind=kind_phys), allocatable, target, dimension(:,:,:,:) :: sfc_var3, phy_var3 character(len=32), allocatable, dimension(:) :: dust12m_name, emi_name, rrfssd_name real(kind=kind_phys), allocatable, target, dimension(:,:,:,:) :: rrfssd_var @@ -758,34 +758,33 @@ pure subroutine fill_Sfcprop_names(Model,sfc_name2,sfc_name3,nvar_s2m,warm_start nt=nt+1 ; sfc_name2(nt) = 'albdifvis_ice' nt=nt+1 ; sfc_name2(nt) = 'albdirnir_ice' nt=nt+1 ; sfc_name2(nt) = 'albdifnir_ice' -! nt=nt+1 ; sfc_name2(nt) = 'sfalb_ice' endif if(Model%cplwav) then - sfc_name2(nvar_s2m) = 'zorlwav' !zorl from wave component + nt=nt+1 ; sfc_name2(nvar_s2m) = 'zorlwav' !zorl from wave component endif - nt = nvar_s2m ! next variable will be at nvar_s2m - + if (Model%nstf_name(1) > 0) then !--- NSSTM inputs only needed when (nstf_name(1) > 0) .and. (nstf_name(2)) == 0) - nt=nt+1 ; sfc_name2(nt) = 'tref' - nt=nt+1 ; sfc_name2(nt) = 'z_c' - nt=nt+1 ; sfc_name2(nt) = 'c_0' - nt=nt+1 ; sfc_name2(nt) = 'c_d' - nt=nt+1 ; sfc_name2(nt) = 'w_0' - nt=nt+1 ; sfc_name2(nt) = 'w_d' - nt=nt+1 ; sfc_name2(nt) = 'xt' - nt=nt+1 ; sfc_name2(nt) = 'xs' - nt=nt+1 ; sfc_name2(nt) = 'xu' - nt=nt+1 ; sfc_name2(nt) = 'xv' - nt=nt+1 ; sfc_name2(nt) = 'xz' - nt=nt+1 ; sfc_name2(nt) = 'zm' - nt=nt+1 ; sfc_name2(nt) = 'xtts' - nt=nt+1 ; sfc_name2(nt) = 'xzts' - nt=nt+1 ; sfc_name2(nt) = 'd_conv' - nt=nt+1 ; sfc_name2(nt) = 'ifd' - nt=nt+1 ; sfc_name2(nt) = 'dt_cool' - nt=nt+1 ; sfc_name2(nt) = 'qrain' + nt=nt+1 ; sfc_name2(nt) = 'tref' + nt=nt+1 ; sfc_name2(nt) = 'z_c' + nt=nt+1 ; sfc_name2(nt) = 'c_0' + nt=nt+1 ; sfc_name2(nt) = 'c_d' + nt=nt+1 ; sfc_name2(nt) = 'w_0' + nt=nt+1 ; sfc_name2(nt) = 'w_d' + nt=nt+1 ; sfc_name2(nt) = 'xt' + nt=nt+1 ; sfc_name2(nt) = 'xs' + nt=nt+1 ; sfc_name2(nt) = 'xu' + nt=nt+1 ; sfc_name2(nt) = 'xv' + nt=nt+1 ; sfc_name2(nt) = 'xz' + nt=nt+1 ; sfc_name2(nt) = 'zm' + nt=nt+1 ; sfc_name2(nt) = 'xtts' + nt=nt+1 ; sfc_name2(nt) = 'xzts' + nt=nt+1 ; sfc_name2(nt) = 'd_conv' + nt=nt+1 ; sfc_name2(nt) = 'ifd' + nt=nt+1 ; sfc_name2(nt) = 'dt_cool' + nt=nt+1 ; sfc_name2(nt) = 'qrain' + endif ! ! Only needed when Noah MP LSM is used - 29 2D ! @@ -862,11 +861,12 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta logical, intent(in) :: warm_start logical, intent(in) :: ignore_rst_cksum !--- local variables - integer :: i, j, k, ix, lsoil, num, nb, i_start, j_start, i_end, j_end, nt + integer :: i, j, k, ix, lsoil, num, nb, i_start, j_start, i_end, j_end, nt, n integer :: isc, iec, jsc, jec, npz, nx, ny integer :: id_restart integer :: nvar_o2, nvar_s2m, nvar_s2o, nvar_s3 integer :: nvar_oro_ls_ss + integer :: nvar_vegfr, nvar_soilfr integer :: nvar_s2r, nvar_s2mp, nvar_s3mp, isnow integer :: nvar_emi, nvar_dust12m, nvar_rrfssd integer, allocatable :: ii1(:), jj1(:) @@ -875,6 +875,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta real(kind=kind_phys), pointer, dimension(:,:,:) :: var3_p1 => NULL() real(kind=kind_phys), pointer, dimension(:,:,:) :: var3_p2 => NULL() real(kind=kind_phys), pointer, dimension(:,:,:) :: var3_p3 => NULL() + real(kind=kind_phys), pointer, dimension(:,:,:) :: var3_fr => NULL() !--- local variables for sncovr calculation integer :: vegtyp logical :: mand @@ -890,7 +891,15 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta nvar_o2 = 19 nvar_oro_ls_ss = 10 - nvar_s2o = 18 + + nvar_vegfr = Model%nvegcat + nvar_soilfr = Model%nsoilcat + + if (Model%nstf_name(1) > 0) then + nvar_s2o = 18 + else + nvar_s2o = 0 + endif if(Model%rrfs_sd) then nvar_dust12m = 5 nvar_rrfssd = 3 @@ -944,42 +953,54 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta !--- allocate the various containers needed for orography data allocate(oro_name2(nvar_o2)) allocate(oro_var2(nx,ny,nvar_o2)) + + allocate(oro_var3v(nx,ny,nvar_vegfr)) + allocate(oro_var3s(nx,ny,nvar_soilfr)) + oro_var2 = -9999._kind_phys - oro_name2(1) = 'stddev' ! hprime(ix,1) - oro_name2(2) = 'convexity' ! hprime(ix,2) - oro_name2(3) = 'oa1' ! hprime(ix,3) - oro_name2(4) = 'oa2' ! hprime(ix,4) - oro_name2(5) = 'oa3' ! hprime(ix,5) - oro_name2(6) = 'oa4' ! hprime(ix,6) - oro_name2(7) = 'ol1' ! hprime(ix,7) - oro_name2(8) = 'ol2' ! hprime(ix,8) - oro_name2(9) = 'ol3' ! hprime(ix,9) - oro_name2(10) = 'ol4' ! hprime(ix,10) - oro_name2(11) = 'theta' ! hprime(ix,11) - oro_name2(12) = 'gamma' ! hprime(ix,12) - oro_name2(13) = 'sigma' ! hprime(ix,13) - oro_name2(14) = 'elvmax' ! hprime(ix,14) - oro_name2(15) = 'orog_filt' ! oro - oro_name2(16) = 'orog_raw' ! oro_uf - oro_name2(17) = 'land_frac' ! land fraction [0:1] + num = 1 ; oro_name2(num) = 'stddev' ! hprime(ix,1) + num = num + 1 ; oro_name2(num) = 'convexity' ! hprime(ix,2) + num = num + 1 ; oro_name2(num) = 'oa1' ! hprime(ix,3) + num = num + 1 ; oro_name2(num) = 'oa2' ! hprime(ix,4) + num = num + 1 ; oro_name2(num) = 'oa3' ! hprime(ix,5) + num = num + 1 ; oro_name2(num) = 'oa4' ! hprime(ix,6) + num = num + 1 ; oro_name2(num) = 'ol1' ! hprime(ix,7) + num = num + 1 ; oro_name2(num) = 'ol2' ! hprime(ix,8) + num = num + 1 ; oro_name2(num) = 'ol3' ! hprime(ix,9) + num = num + 1 ; oro_name2(num) = 'ol4' ! hprime(ix,10) + num = num + 1 ; oro_name2(num) = 'theta' ! hprime(ix,11) + num = num + 1 ; oro_name2(num) = 'gamma' ! hprime(ix,12) + num = num + 1 ; oro_name2(num) = 'sigma' ! hprime(ix,13) + num = num + 1 ; oro_name2(num) = 'elvmax' ! hprime(ix,14) + num = num + 1 ; oro_name2(num) = 'orog_filt' ! oro + num = num + 1 ; oro_name2(num) = 'orog_raw' ! oro_uf + num = num + 1 ; oro_name2(num) = 'land_frac' ! land fraction [0:1] !--- variables below here are optional - oro_name2(18) = 'lake_frac' ! lake fraction [0:1] - oro_name2(19) = 'lake_depth' ! lake depth(m) + num = num + 1 ; oro_name2(num) = 'lake_frac' ! lake fraction [0:1] + num = num + 1 ; oro_name2(num) = 'lake_depth' ! lake depth(m) !--- register axis call register_axis( Oro_restart, "lon", 'X' ) call register_axis( Oro_restart, "lat", 'Y' ) !--- register the 2D fields - do num = 1,nvar_o2 - var2_p => oro_var2(:,:,num) - if (trim(oro_name2(num)) == 'lake_frac' .or. trim(oro_name2(num)) == 'lake_depth') then - call register_restart_field(Oro_restart, oro_name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.) + do n = 1,num + var2_p => oro_var2(:,:,n) + if (trim(oro_name2(n)) == 'lake_frac' .or. trim(oro_name2(n)) == 'lake_depth' ) then + call register_restart_field(Oro_restart, oro_name2(n), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.) else - call register_restart_field(Oro_restart, oro_name2(num), var2_p, dimensions=(/'lat','lon'/)) + call register_restart_field(Oro_restart, oro_name2(n), var2_p, dimensions=(/'lat','lon'/)) endif enddo nullify(var2_p) + + !--- register 3D vegetation and soil fractions + var3_fr => oro_var3v(:,:,:) + call register_restart_field(Oro_restart, 'vegetation_type_pct', var3_fr, dimensions=(/'num_veg_cat','lat ','lon '/) , is_optional=.true.) + var3_fr => oro_var3s(:,:,:) + call register_restart_field(Oro_restart, 'soil_type_pct', var3_fr, dimensions=(/'num_soil_cat','lat ','lon '/) , is_optional=.true.) + nullify(var3_fr) + endif !--- read the orography restart/data @@ -990,7 +1011,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta !--- copy data into GFS containers -!$omp parallel do default(shared) private(i, j, nb, ix) +!$omp parallel do default(shared) private(i, j, nb, ix, num) do nb = 1, Atm_block%nblks !--- 2D variables do ix = 1, Atm_block%blksz(nb) @@ -999,32 +1020,43 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta !--- stddev ! Sfcprop(nb)%hprim(ix) = oro_var2(i,j,1) !--- hprime(1:14) - Sfcprop(nb)%hprime(ix,1) = oro_var2(i,j,1) - Sfcprop(nb)%hprime(ix,2) = oro_var2(i,j,2) - Sfcprop(nb)%hprime(ix,3) = oro_var2(i,j,3) - Sfcprop(nb)%hprime(ix,4) = oro_var2(i,j,4) - Sfcprop(nb)%hprime(ix,5) = oro_var2(i,j,5) - Sfcprop(nb)%hprime(ix,6) = oro_var2(i,j,6) - Sfcprop(nb)%hprime(ix,7) = oro_var2(i,j,7) - Sfcprop(nb)%hprime(ix,8) = oro_var2(i,j,8) - Sfcprop(nb)%hprime(ix,9) = oro_var2(i,j,9) - Sfcprop(nb)%hprime(ix,10) = oro_var2(i,j,10) - Sfcprop(nb)%hprime(ix,11) = oro_var2(i,j,11) - Sfcprop(nb)%hprime(ix,12) = oro_var2(i,j,12) - Sfcprop(nb)%hprime(ix,13) = oro_var2(i,j,13) - Sfcprop(nb)%hprime(ix,14) = oro_var2(i,j,14) + num = 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%hprime(ix,num) = oro_var2(i,j,num) !--- oro - Sfcprop(nb)%oro(ix) = oro_var2(i,j,15) - !--- oro_uf - Sfcprop(nb)%oro_uf(ix) = oro_var2(i,j,16) + num = num + 1 ; Sfcprop(nb)%oro(ix) = oro_var2(i,j,num) + num = num + 1 ; Sfcprop(nb)%oro_uf(ix) = oro_var2(i,j,num) Sfcprop(nb)%landfrac(ix) = -9999.0 Sfcprop(nb)%lakefrac(ix) = -9999.0 - Sfcprop(nb)%landfrac(ix) = oro_var2(i,j,17) !land frac [0:1] - Sfcprop(nb)%lakefrac(ix) = oro_var2(i,j,18) !lake frac [0:1] + num = num + 1 ; Sfcprop(nb)%landfrac(ix) = oro_var2(i,j,num) !land frac [0:1] + num = num + 1 ; Sfcprop(nb)%lakefrac(ix) = oro_var2(i,j,num) !lake frac [0:1] + num = num + 1 ; Sfcprop(nb)%lakedepth(ix) = oro_var2(i,j,num) !lake depth [m] !YWu + + Sfcprop(nb)%vegtype_frac(ix,:) = -9999.0 + Sfcprop(nb)%soiltype_frac(ix,:) = -9999.0 + + Sfcprop(nb)%vegtype_frac(ix,:) = oro_var3v(i,j,:) ! vegetation type fractions, [0:1] + Sfcprop(nb)%soiltype_frac(ix,:) = oro_var3s(i,j,:) ! soil type fractions, [0:1] - Sfcprop(nb)%lakedepth(ix) = oro_var2(i,j,19) !lake depth [m] !YWu + !do n=1,nvar_vegfr + ! if (Sfcprop(nb)%vegtype_frac(ix,n) > 0.) print *,'Sfcprop(nb)%vegtype_frac(ix,n)',Sfcprop(nb)%vegtype_frac(ix,n),n + !enddo + !do n=1,nvar_soilfr + ! if (Sfcprop(nb)%soiltype_frac(ix,n) > 0.) print *,'Sfcprop(nb)%soiltype_frac(ix,n)',Sfcprop(nb)%soiltype_frac(ix,n),n + !enddo enddo enddo @@ -1040,6 +1072,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta !--- deallocate containers and free restart container deallocate(oro_name2, oro_var2) + deallocate(oro_var3v) + deallocate(oro_var3s) if_smoke: if(Model%rrfs_sd) then ! for RRFS-SD @@ -1590,7 +1624,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta ! call copy_to_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sfalb_ice) endif if(Model%cplwav) then - nt = nvar_s2m-1 ! Next item will be at nvar_s2m + !tgs - the following line is a bug. It should be nt = nt + !nt = nvar_s2m-1 ! Next item will be at nvar_s2m call copy_to_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorlwav) !--- (zorl from wave model) else Sfcprop(nb)%zorlwav = Sfcprop(nb)%zorlw @@ -1703,7 +1738,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta ! !--- NSSTM variables - nt = nvar_s2m + !tgs - the following line is a bug that will show if(Model%cplwav) = true + !nt = nvar_s2m if (Model%nstf_name(1) > 0) then if (Model%nstf_name(2) == 1) then ! nsst spinup !--- nsstm tref @@ -1746,8 +1782,6 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta call copy_to_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%dt_cool) !--- nsstm dt_cool call copy_to_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%qrain) !--- nsstm qrain endif - else - nt = nt + 18 endif if (Model%lsm == Model%lsm_ruc .and. warm_start) then @@ -2134,7 +2168,11 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta ! nvar2m = nvar2m + 5 endif if (Model%cplwav) nvar2m = nvar2m + 1 - nvar2o = 18 + if (Model%nstf_name(1) > 0) then + nvar2o = 18 + else + nvar2o = 0 + endif if (Model%lsm == Model%lsm_ruc) then if (Model%rdlai) then nvar2r = 13 @@ -2499,7 +2537,7 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sfalb_ice) if (Model%rdlai) then call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xlaixy) - endif + endif else if (Model%lsm == Model%lsm_noahmp) then !--- Extra Noah MP variables call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snowxy) diff --git a/io/FV3GFS_restart_io.F90 b/io/FV3GFS_restart_io.F90 new file mode 100644 index 000000000..da641c04e --- /dev/null +++ b/io/FV3GFS_restart_io.F90 @@ -0,0 +1,920 @@ +module FV3GFS_restart_io_mod + + use esmf + use block_control_mod, only: block_control_type + use GFS_typedefs, only: GFS_sfcprop_type, GFS_control_type, kind_phys + use GFS_restart, only: GFS_restart_type + + implicit none + private + + real(kind_phys), parameter:: zero = 0.0, one = 1.0 + integer, parameter :: r8 = kind_phys + + integer :: nvar2d, nvar3d, npz + real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: phy_var2 + real(kind=kind_phys), allocatable, target, dimension(:,:,:,:) :: phy_var3 + character(len=32),dimension(:),allocatable :: phy_var2_names, phy_var3_names + + integer :: nvar2m, nvar2o, nvar3, nvar2r, nvar2mp, nvar3mp + real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: sfc_var2 + real(kind=kind_phys), allocatable, target, dimension(:,:,:) :: sfc_var3ice + real(kind=kind_phys), allocatable, target, dimension(:,:,:,:) :: sfc_var3, sfc_var3sn,sfc_var3eq,sfc_var3zn + character(len=32),allocatable,dimension(:) :: sfc_name2, sfc_name3 + + type(ESMF_FieldBundle) :: phy_bundle, sfc_bundle + + public FV3GFS_restart_register + + public fv_phy_restart_output + public fv_phy_restart_bundle_setup + + public fv_sfc_restart_output + public fv_sfc_restart_bundle_setup + + interface copy_from_GFS_Data + module procedure copy_from_GFS_Data_2d_phys2phys, & + copy_from_GFS_Data_3d_phys2phys, & + copy_from_GFS_Data_2d_int2phys, & + copy_from_GFS_Data_3d_int2phys, & + copy_from_GFS_Data_2d_stack_phys2phys + end interface + + contains + + subroutine FV3GFS_restart_register (Sfcprop, GFS_restart, Atm_block, Model) + + ! this subroutine must allocate all data buffers and set the variable names + ! for both 'phy' and 'sfc' restart bundles + + implicit none + + type(GFS_sfcprop_type), intent(in) :: Sfcprop(:) + type(GFS_restart_type), intent(in) :: GFS_Restart + type(block_control_type), intent(in) :: Atm_block + type(GFS_control_type), intent(in) :: Model + + integer :: isc, iec, jsc, jec, nx, ny + integer :: num + + isc = Atm_block%isc + iec = Atm_block%iec + jsc = Atm_block%jsc + jec = Atm_block%jec + npz = Atm_block%npz + nx = (iec - isc + 1) + ny = (jec - jsc + 1) + + !--------------- phy + nvar2d = GFS_Restart%num2d + nvar3d = GFS_Restart%num3d + + allocate (phy_var2(nx,ny,nvar2d), phy_var2_names(nvar2d)) + allocate (phy_var3(nx,ny,npz,nvar3d), phy_var3_names(nvar3d)) + phy_var2 = zero + phy_var3 = zero + do num = 1,nvar2d + phy_var2_names(num) = trim(GFS_Restart%name2d(num)) + enddo + do num = 1,nvar3d + phy_var3_names(num) = trim(GFS_Restart%name3d(num)) + enddo + + !--------------- sfc + nvar2m = 48 + if (Model%use_cice_alb .or. Model%lsm == Model%lsm_ruc) then + nvar2m = nvar2m + 4 +! nvar2m = nvar2m + 5 + endif + if (Model%cplwav) nvar2m = nvar2m + 1 + nvar2o = 18 + if (Model%lsm == Model%lsm_ruc) then + if (Model%rdlai) then + nvar2r = 13 + else + nvar2r = 12 + endif + nvar3 = 5 + else + nvar2r = 0 + nvar3 = 3 + endif + nvar2mp = 0 + nvar3mp = 0 + if (Model%lsm == Model%lsm_noahmp) then + nvar2mp = 29 + nvar3mp = 5 + endif + + !--- allocate the various containers needed for restarts + allocate(sfc_name2(nvar2m+nvar2o+nvar2mp+nvar2r)) + allocate(sfc_var2(nx,ny,nvar2m+nvar2o+nvar2mp+nvar2r)) + allocate(sfc_name3(0:nvar3+nvar3mp)) + allocate(sfc_var3ice(nx,ny,Model%kice)) + + if (Model%lsm == Model%lsm_noah .or. Model%lsm == Model%lsm_noahmp) then + allocate(sfc_var3(nx,ny,Model%lsoil,nvar3)) + elseif (Model%lsm == Model%lsm_ruc) then + allocate(sfc_var3(nx,ny,Model%lsoil_lsm,nvar3)) + endif + + sfc_var2 = -9999.0_r8 + sfc_var3 = -9999.0_r8 + sfc_var3ice= -9999.0_r8 + + if (Model%lsm == Model%lsm_noahmp) then + allocate(sfc_var3sn(nx,ny,-2:0,4:6)) + allocate(sfc_var3eq(nx,ny,1:4,7:7)) + allocate(sfc_var3zn(nx,ny,-2:4,8:8)) + + sfc_var3sn = -9999.0_r8 + sfc_var3eq = -9999.0_r8 + sfc_var3zn = -9999.0_r8 + endif + + call fill_Sfcprop_names(Model,sfc_name2,sfc_name3,nvar2m,.true.) + + sfc_name3(0) = 'tiice' + + if (Model%lsm == Model%lsm_noah .or. Model%lsm == Model%lsm_noahmp) then + sfc_name3(1) = 'stc' + sfc_name3(2) = 'smc' + sfc_name3(3) = 'slc' + if (Model%lsm == Model%lsm_noahmp) then + sfc_name3(4) = 'snicexy' + sfc_name3(5) = 'snliqxy' + sfc_name3(6) = 'tsnoxy' + sfc_name3(7) = 'smoiseq' + sfc_name3(8) = 'zsnsoxy' + endif + else if (Model%lsm == Model%lsm_ruc) then + sfc_name3(1) = 'tslb' + sfc_name3(2) = 'smois' + sfc_name3(3) = 'sh2o' + sfc_name3(4) = 'smfr' + sfc_name3(5) = 'flfr' + end if + + end subroutine FV3GFS_restart_register + + subroutine fv_phy_restart_output(GFS_Restart, Atm_block) + + implicit none + + type(GFS_restart_type), intent(in) :: GFS_Restart + type(block_control_type), intent(in) :: Atm_block + +!*** local variables + integer :: i, j, k, n + integer :: nb, ix, num + integer :: isc, iec, jsc, jec, npz, nx, ny + integer(8) :: rchk + + isc = Atm_block%isc + iec = Atm_block%iec + jsc = Atm_block%jsc + jec = Atm_block%jec + npz = Atm_block%npz + nx = (iec - isc + 1) + ny = (jec - jsc + 1) + + !--- register the restart fields + if (.not. allocated(phy_var2)) then + write(0,*)'phy_var2 must be allocated' + endif + if (.not. allocated(phy_var3)) then + write(0,*)'phy_var3 must be allocated' + endif + + !--- 2D variables + do num = 1,nvar2d + do nb = 1,Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + i = Atm_block%index(nb)%ii(ix) - isc + 1 + j = Atm_block%index(nb)%jj(ix) - jsc + 1 + phy_var2(i,j,num) = GFS_Restart%data(nb,num)%var2p(ix) + enddo + enddo + enddo + + !--- 3D variables + do num = 1,nvar3d + do nb = 1,Atm_block%nblks + do k=1,npz + do ix = 1, Atm_block%blksz(nb) + i = Atm_block%index(nb)%ii(ix) - isc + 1 + j = Atm_block%index(nb)%jj(ix) - jsc + 1 + phy_var3(i,j,k,num) = GFS_Restart%data(nb,num)%var3p(ix,k) + enddo + enddo + enddo + enddo + + end subroutine fv_phy_restart_output + + subroutine fv_sfc_restart_output(Sfcprop, Atm_block, Model) + !--- interface variable definitions + implicit none + + type(GFS_sfcprop_type), intent(in) :: Sfcprop(:) + type(block_control_type), intent(in) :: Atm_block + type(GFS_control_type), intent(in) :: Model + + integer :: i, j, k, nb, ix, lsoil, num, nt + integer :: isc, iec, jsc, jec, npz, nx, ny + integer, allocatable :: ii1(:), jj1(:) + real(kind_phys) :: ice + integer :: is, ie + + isc = Atm_block%isc + iec = Atm_block%iec + jsc = Atm_block%jsc + jec = Atm_block%jec + npz = Atm_block%npz + nx = (iec - isc + 1) + ny = (jec - jsc + 1) + +!$omp parallel do default(shared) private(i, j, nb, ix, nt, ii1, jj1, lsoil, k, ice) + block_loop: do nb = 1, Atm_block%nblks + allocate(ii1(Atm_block%blksz(nb))) + allocate(jj1(Atm_block%blksz(nb))) + ii1=Atm_block%index(nb)%ii - isc + 1 + jj1=Atm_block%index(nb)%jj - jsc + 1 + + nt=0 + + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%slmsk) !--- slmsk + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tsfco) !--- tsfc (tsea in sfc file) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%weasd) !--- weasd (sheleg in sfc file) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tg3) !--- tg3 + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorl) !--- zorl + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%alvsf) !--- alvsf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%alvwf) !--- alvwf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%alnsf) !--- alnsf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%alnwf) !--- alnwf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%facsf) !--- facsf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%facwf) !--- facwf + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%vfrac) !--- vfrac + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%canopy)!--- canopy + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%f10m) !--- f10m + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%t2m) !--- t2m + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%q2m) !--- q2m + + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%vtype) !--- vtype + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%stype) !--- stype + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%uustar)!--- uustar + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%ffmm) !--- ffmm + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%ffhh) !--- ffhh + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%hice) !--- hice + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%fice) !--- fice + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tisfc) !--- tisfc + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tprcp) !--- tprcp + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%srflag)!--- srflag + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snowd) !--- snowd (snwdph in the file) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%shdmin)!--- shdmin + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%shdmax)!--- shdmax + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%slope) !--- slope + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snoalb)!--- snoalb + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sncovr) !--- sncovr + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snodl) !--- snodl (snowd on land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%weasdl) !--- weasdl (weasd on land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tsfc) !--- tsfc composite + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tsfcl) !--- tsfcl (temp on land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorlw) !--- zorl (zorl on water) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorll) !--- zorll (zorl on land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorli) !--- zorli (zorl on ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdirvis_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdirnir_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdifvis_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdifnir_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%emis_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%emis_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sncovr_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snodi) !--- snodi (snowd on ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%weasdi) !--- weasdi (weasd on ice) + if (Model%use_cice_alb .or. Model%lsm == Model%lsm_ruc) then + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdirvis_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdifvis_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdirnir_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%albdifnir_ice) +! sfc_var2(i,j,53) = Sfcprop(nb)%sfalb_ice(ix) + endif + if (Model%cplwav) then + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zorlwav) !--- zorlwav (zorl from wav) + endif + !--- NSSTM variables + if (Model%nstf_name(1) > 0) then + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tref) !--- nsstm tref + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%z_c) !--- nsstm z_c + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%c_0) !--- nsstm c_0 + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%c_d) !--- nsstm c_d + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%w_0) !--- nsstm w_0 + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%w_d) !--- nsstm w_d + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xt) !--- nsstm xt + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xs) !--- nsstm xs + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xu) !--- nsstm xu + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xv) !--- nsstm xv + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xz) !--- nsstm xz + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zm) !--- nsstm zm + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xtts) !--- nsstm xtts + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xzts) !--- nsstm xzts + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%d_conv) !--- nsstm d_conv + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%ifd) !--- nsstm ifd + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%dt_cool)!--- nsstm dt_cool + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%qrain) !--- nsstm qrain + + ! FIXME convert negative zero (-0.0) to zero (0.0) + do j=1,ny + do i=1,nx + if(sfc_var2(i,j,nt) == 0.0) sfc_var2(i,j,nt) = 0.0 + end do + end do + endif + + if (Model%lsm == Model%lsm_ruc) then + !--- Extra RUC variables + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%wetness) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%clw_surf_land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%clw_surf_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%qwv_surf_land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%qwv_surf_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tsnow_land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tsnow_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snowfallac_land) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snowfallac_ice) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sfalb_lnd) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sfalb_lnd_bck) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sfalb_ice) + if (Model%rdlai) then + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xlaixy) + endif + else if (Model%lsm == Model%lsm_noahmp) then + !--- Extra Noah MP variables + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%snowxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tvxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tgxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%canicexy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%canliqxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%eahxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%tahxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%cmxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%chxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%fwetxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%sneqvoxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%alboldxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%qsnowxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%wslakexy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%zwtxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%waxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%wtxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%lfmassxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%rtmassxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%stmassxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%woodxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%stblcpxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%fastcpxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xsaixy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%xlaixy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%taussxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%smcwtdxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%deeprechxy) + call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,sfc_var2,Sfcprop(nb)%rechxy) + endif + + do k = 1,Model%kice + do ix = 1, Atm_block%blksz(nb) + ice=Sfcprop(nb)%tiice(ix,k) + if(ice phy_var2(:,:,num) + call create_2d_field_and_add_to_bundle(temp_r2d, trim(phy_var2_names(num)), trim(outputfile), grid, bundle) + enddo + + do num = 1,nvar3d + temp_r3d => phy_var3(:,:,:,num) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(phy_var3_names(num)), "zaxis_1", npz, trim(outputfile), grid, bundle) + enddo + + end subroutine fv_phy_restart_bundle_setup + + subroutine fv_sfc_restart_bundle_setup(bundle, grid, Model, rc) +! +!------------------------------------------------------------- +!*** set esmf bundle for sfc restart fields +!------------------------------------------------------------ +! + use esmf + + implicit none + + type(ESMF_FieldBundle),intent(inout) :: bundle + type(ESMF_Grid),intent(inout) :: grid + type(GFS_control_type), intent(in) :: Model + integer,intent(out) :: rc + +!*** local variables + integer i, j, k, n + character(128) :: sfcbdl_name + type(ESMF_Field) :: field + character(128) :: outputfile + real(kind_phys),dimension(:,:),pointer :: temp_r2d + real(kind_phys),dimension(:,:,:),pointer :: temp_r3d + + integer :: num + + if (.not. allocated(sfc_var2)) then + write(0,*)'ERROR sfc_var2, NOT allocated' + endif + if (.not. allocated(sfc_var3)) then + write(0,*)'ERROR sfc_var3 NOT allocated' + endif + + sfc_bundle = bundle + + call ESMF_FieldBundleGet(bundle, name=sfcbdl_name,rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + outputfile = trim(sfcbdl_name) + +!*** add esmf fields + + do num = 1,nvar2m + temp_r2d => sfc_var2(:,:,num) + call create_2d_field_and_add_to_bundle(temp_r2d, trim(sfc_name2(num)), outputfile, grid, bundle) + enddo + + if (Model%nstf_name(1) > 0) then + do num = nvar2m+1,nvar2m+nvar2o + temp_r2d => sfc_var2(:,:,num) + call create_2d_field_and_add_to_bundle(temp_r2d, trim(sfc_name2(num)), outputfile, grid, bundle) + enddo + endif + + if (Model%lsm == Model%lsm_ruc) then ! nvar2mp =0 + do num = nvar2m+nvar2o+1, nvar2m+nvar2o+nvar2r + temp_r2d => sfc_var2(:,:,num) + call create_2d_field_and_add_to_bundle(temp_r2d, trim(sfc_name2(num)), outputfile, grid, bundle) + enddo + else if (Model%lsm == Model%lsm_noahmp) then ! nvar2r =0 + do num = nvar2m+nvar2o+1,nvar2m+nvar2o+nvar2mp + temp_r2d => sfc_var2(:,:,num) + call create_2d_field_and_add_to_bundle(temp_r2d, trim(sfc_name2(num)), outputfile, grid, bundle) + enddo + endif + + temp_r3d => sfc_var3ice(:,:,:) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(0)), "zaxis_1", Model%kice, trim(outputfile), grid, bundle) + + if(Model%lsm == Model%lsm_ruc) then + do num = 1,nvar3 + temp_r3d => sfc_var3(:,:,:,num) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(num)), "zaxis_1", Model%kice, trim(outputfile), grid, bundle) + enddo + else + do num = 1,nvar3 + temp_r3d => sfc_var3(:,:,:,num) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(num)), "zaxis_2", Model%lsoil, trim(outputfile), grid, bundle) + enddo + endif + + if (Model%lsm == Model%lsm_noahmp) then + do num = nvar3+1,nvar3+3 + temp_r3d => sfc_var3sn(:,:,:,num) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(num)), "zaxis_3", 3, trim(outputfile), grid, bundle) + enddo + + temp_r3d => sfc_var3eq(:,:,:,7) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(7)), "zaxis_2", Model%lsoil, trim(outputfile), grid, bundle) + + temp_r3d => sfc_var3zn(:,:,:,8) + call create_3d_field_and_add_to_bundle(temp_r3d, trim(sfc_name3(8)), "zaxis_4", 7, trim(outputfile), grid, bundle) + endif ! lsm = lsm_noahmp + + end subroutine fv_sfc_restart_bundle_setup + + subroutine create_2d_field_and_add_to_bundle(temp_r2d, field_name, outputfile, grid, bundle) + + use esmf + + implicit none + + real(kind_phys), dimension(:,:), pointer, intent(in) :: temp_r2d + character(len=*), intent(in) :: field_name + character(len=*), intent(in) :: outputfile + type(ESMF_Grid), intent(in) :: grid + type(ESMF_FieldBundle), intent(inout) :: bundle + + type(ESMF_Field) :: field + + integer :: rc, i + + field = ESMF_FieldCreate(grid, temp_r2d, datacopyflag=ESMF_DATACOPY_REFERENCE, & + name=trim(field_name), indexFlag=ESMF_INDEX_DELOCAL, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__, file=__FILE__)) & + call ESMF_Finalize(endflag=ESMF_END_ABORT) + + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", name='output_file', value=trim(outputfile), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) + + call ESMF_FieldBundleAdd(bundle, (/field/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) + + end subroutine create_2d_field_and_add_to_bundle + + subroutine create_3d_field_and_add_to_bundle(temp_r3d, field_name, axis_name, num_levels, outputfile, grid, bundle) + + use esmf + + implicit none + + real(kind_phys), dimension(:,:,:), pointer, intent(in) :: temp_r3d + character(len=*), intent(in) :: field_name + character(len=*), intent(in) :: axis_name + integer, intent(in) :: num_levels + character(len=*), intent(in) :: outputfile + type(ESMF_Grid), intent(in) :: grid + type(ESMF_FieldBundle), intent(inout) :: bundle + + type(ESMF_Field) :: field + + integer :: rc, i + + field = ESMF_FieldCreate(grid, temp_r3d, datacopyflag=ESMF_DATACOPY_REFERENCE, & + name=trim(field_name), indexFlag=ESMF_INDEX_DELOCAL, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__, file=__FILE__)) & + call ESMF_Finalize(endflag=ESMF_END_ABORT) + + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", name='output_file', value=trim(outputfile), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) + + call add_zaxis_to_field(field, axis_name, num_levels) + + call ESMF_FieldBundleAdd(bundle, (/field/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) + + end subroutine create_3d_field_and_add_to_bundle + + subroutine add_zaxis_to_field(field, axis_name, num_levels) + + use esmf + + implicit none + + type(ESMF_Field), intent(inout) :: field + character(len=*), intent(in) :: axis_name + integer, intent(in) :: num_levels + + real(kind_phys), allocatable, dimension(:) :: buffer + integer :: rc, i + + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & + name="ESMF:ungridded_dim_labels", valueList=(/trim(axis_name)/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + allocate( buffer(num_levels) ) + do i=1, num_levels + buffer(i)=i + end do + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3-dim", & + name=trim(axis_name), valueList=buffer, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + deallocate(buffer) + + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3-dim", & + name=trim(axis_name)//"cartesian_axis", value="Z", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + end subroutine add_zaxis_to_field + + pure subroutine fill_Sfcprop_names(Model,sfc_name2,sfc_name3,nvar_s2m,warm_start) + implicit none + type(GFS_control_type), intent(in) :: Model + integer, intent(in) :: nvar_s2m + character(len=32),intent(out) :: sfc_name2(:), sfc_name3(:) + logical, intent(in) :: warm_start + integer :: nt + + !--- names of the 2D variables to save + nt=0 + nt=nt+1 ; sfc_name2(nt) = 'slmsk' + nt=nt+1 ; sfc_name2(nt) = 'tsea' !tsfc + nt=nt+1 ; sfc_name2(nt) = 'sheleg' !weasd + nt=nt+1 ; sfc_name2(nt) = 'tg3' + nt=nt+1 ; sfc_name2(nt) = 'zorl' + nt=nt+1 ; sfc_name2(nt) = 'alvsf' + nt=nt+1 ; sfc_name2(nt) = 'alvwf' + nt=nt+1 ; sfc_name2(nt) = 'alnsf' + nt=nt+1 ; sfc_name2(nt) = 'alnwf' + nt=nt+1 ; sfc_name2(nt) = 'facsf' + nt=nt+1 ; sfc_name2(nt) = 'facwf' + nt=nt+1 ; sfc_name2(nt) = 'vfrac' + nt=nt+1 ; sfc_name2(nt) = 'canopy' + nt=nt+1 ; sfc_name2(nt) = 'f10m' + nt=nt+1 ; sfc_name2(nt) = 't2m' + nt=nt+1 ; sfc_name2(nt) = 'q2m' + nt=nt+1 ; sfc_name2(nt) = 'vtype' + nt=nt+1 ; sfc_name2(nt) = 'stype' + nt=nt+1 ; sfc_name2(nt) = 'uustar' + nt=nt+1 ; sfc_name2(nt) = 'ffmm' + nt=nt+1 ; sfc_name2(nt) = 'ffhh' + nt=nt+1 ; sfc_name2(nt) = 'hice' + nt=nt+1 ; sfc_name2(nt) = 'fice' + nt=nt+1 ; sfc_name2(nt) = 'tisfc' + nt=nt+1 ; sfc_name2(nt) = 'tprcp' + nt=nt+1 ; sfc_name2(nt) = 'srflag' + nt=nt+1 ; sfc_name2(nt) = 'snwdph' !snowd + nt=nt+1 ; sfc_name2(nt) = 'shdmin' + nt=nt+1 ; sfc_name2(nt) = 'shdmax' + nt=nt+1 ; sfc_name2(nt) = 'slope' + nt=nt+1 ; sfc_name2(nt) = 'snoalb' + !--- variables below here are optional + nt=nt+1 ; sfc_name2(nt) = 'sncovr' + nt=nt+1 ; sfc_name2(nt) = 'snodl' !snowd on land portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'weasdl'!weasd on land portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'tsfc' !tsfc composite + nt=nt+1 ; sfc_name2(nt) = 'tsfcl' !temp on land portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'zorlw' !zorl on water portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'zorll' !zorl on land portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'zorli' !zorl on ice portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'albdirvis_lnd' + nt=nt+1 ; sfc_name2(nt) = 'albdirnir_lnd' + nt=nt+1 ; sfc_name2(nt) = 'albdifvis_lnd' + nt=nt+1 ; sfc_name2(nt) = 'albdifnir_lnd' + nt=nt+1 ; sfc_name2(nt) = 'emis_lnd' + nt=nt+1 ; sfc_name2(nt) = 'emis_ice' + nt=nt+1 ; sfc_name2(nt) = 'sncovr_ice' + nt=nt+1 ; sfc_name2(nt) = 'snodi' ! snowd on ice portion of a cell + nt=nt+1 ; sfc_name2(nt) = 'weasdi'! weasd on ice portion of a cell + + if (Model%use_cice_alb .or. Model%lsm == Model%lsm_ruc) then + nt=nt+1 ; sfc_name2(nt) = 'albdirvis_ice' + nt=nt+1 ; sfc_name2(nt) = 'albdifvis_ice' + nt=nt+1 ; sfc_name2(nt) = 'albdirnir_ice' + nt=nt+1 ; sfc_name2(nt) = 'albdifnir_ice' +! nt=nt+1 ; sfc_name2(nt) = 'sfalb_ice' + endif + + if(Model%cplwav) then + sfc_name2(nvar_s2m) = 'zorlwav' !zorl from wave component + endif + + nt = nvar_s2m ! next variable will be at nvar_s2m + + !--- NSSTM inputs only needed when (nstf_name(1) > 0) .and. (nstf_name(2)) == 0) + nt=nt+1 ; sfc_name2(nt) = 'tref' + nt=nt+1 ; sfc_name2(nt) = 'z_c' + nt=nt+1 ; sfc_name2(nt) = 'c_0' + nt=nt+1 ; sfc_name2(nt) = 'c_d' + nt=nt+1 ; sfc_name2(nt) = 'w_0' + nt=nt+1 ; sfc_name2(nt) = 'w_d' + nt=nt+1 ; sfc_name2(nt) = 'xt' + nt=nt+1 ; sfc_name2(nt) = 'xs' + nt=nt+1 ; sfc_name2(nt) = 'xu' + nt=nt+1 ; sfc_name2(nt) = 'xv' + nt=nt+1 ; sfc_name2(nt) = 'xz' + nt=nt+1 ; sfc_name2(nt) = 'zm' + nt=nt+1 ; sfc_name2(nt) = 'xtts' + nt=nt+1 ; sfc_name2(nt) = 'xzts' + nt=nt+1 ; sfc_name2(nt) = 'd_conv' + nt=nt+1 ; sfc_name2(nt) = 'ifd' + nt=nt+1 ; sfc_name2(nt) = 'dt_cool' + nt=nt+1 ; sfc_name2(nt) = 'qrain' +! +! Only needed when Noah MP LSM is used - 29 2D +! + if (Model%lsm == Model%lsm_noahmp) then + nt=nt+1 ; sfc_name2(nt) = 'snowxy' + nt=nt+1 ; sfc_name2(nt) = 'tvxy' + nt=nt+1 ; sfc_name2(nt) = 'tgxy' + nt=nt+1 ; sfc_name2(nt) = 'canicexy' + nt=nt+1 ; sfc_name2(nt) = 'canliqxy' + nt=nt+1 ; sfc_name2(nt) = 'eahxy' + nt=nt+1 ; sfc_name2(nt) = 'tahxy' + nt=nt+1 ; sfc_name2(nt) = 'cmxy' + nt=nt+1 ; sfc_name2(nt) = 'chxy' + nt=nt+1 ; sfc_name2(nt) = 'fwetxy' + nt=nt+1 ; sfc_name2(nt) = 'sneqvoxy' + nt=nt+1 ; sfc_name2(nt) = 'alboldxy' + nt=nt+1 ; sfc_name2(nt) = 'qsnowxy' + nt=nt+1 ; sfc_name2(nt) = 'wslakexy' + nt=nt+1 ; sfc_name2(nt) = 'zwtxy' + nt=nt+1 ; sfc_name2(nt) = 'waxy' + nt=nt+1 ; sfc_name2(nt) = 'wtxy' + nt=nt+1 ; sfc_name2(nt) = 'lfmassxy' + nt=nt+1 ; sfc_name2(nt) = 'rtmassxy' + nt=nt+1 ; sfc_name2(nt) = 'stmassxy' + nt=nt+1 ; sfc_name2(nt) = 'woodxy' + nt=nt+1 ; sfc_name2(nt) = 'stblcpxy' + nt=nt+1 ; sfc_name2(nt) = 'fastcpxy' + nt=nt+1 ; sfc_name2(nt) = 'xsaixy' + nt=nt+1 ; sfc_name2(nt) = 'xlaixy' + nt=nt+1 ; sfc_name2(nt) = 'taussxy' + nt=nt+1 ; sfc_name2(nt) = 'smcwtdxy' + nt=nt+1 ; sfc_name2(nt) = 'deeprechxy' + nt=nt+1 ; sfc_name2(nt) = 'rechxy' + else if (Model%lsm == Model%lsm_ruc .and. warm_start) then + nt=nt+1 ; sfc_name2(nt) = 'wetness' + nt=nt+1 ; sfc_name2(nt) = 'clw_surf_land' + nt=nt+1 ; sfc_name2(nt) = 'clw_surf_ice' + nt=nt+1 ; sfc_name2(nt) = 'qwv_surf_land' + nt=nt+1 ; sfc_name2(nt) = 'qwv_surf_ice' + nt=nt+1 ; sfc_name2(nt) = 'tsnow_land' + nt=nt+1 ; sfc_name2(nt) = 'tsnow_ice' + nt=nt+1 ; sfc_name2(nt) = 'snowfall_acc_land' + nt=nt+1 ; sfc_name2(nt) = 'snowfall_acc_ice' + nt=nt+1 ; sfc_name2(nt) = 'sfalb_lnd' + nt=nt+1 ; sfc_name2(nt) = 'sfalb_lnd_bck' + nt=nt+1 ; sfc_name2(nt) = 'sfalb_ice' + if (Model%rdlai) then + nt=nt+1 ; sfc_name2(nt) = 'lai' + endif + else if (Model%lsm == Model%lsm_ruc .and. Model%rdlai) then + nt=nt+1 ; sfc_name2(nt) = 'lai' + endif + end subroutine fill_sfcprop_names + + pure subroutine copy_from_GFS_Data_2d_phys2phys(ii1,jj1,isc,jsc,nt,var2d,var_block) + implicit none + integer, intent(in) :: ii1(:), jj1(:), isc, jsc + integer, intent(inout) :: nt + real(kind=kind_phys), intent(in) :: var_block(:) + real(kind=kind_phys), intent(out) :: var2d(:,:,:) + integer ix + + nt=nt+1 + do ix=1,size(var_block) + var2d(ii1(ix),jj1(ix),nt) = var_block(ix) + enddo + end subroutine copy_from_GFS_Data_2d_phys2phys + + pure subroutine copy_from_GFS_Data_3d_phys2phys(ii1,jj1,isc,jsc,nt,var3d,var_block) + implicit none + integer, intent(in) :: ii1(:), jj1(:), isc, jsc + integer, intent(inout) :: nt + real(kind=kind_phys), intent(in) :: var_block(:,:) + real(kind=kind_phys), intent(out) :: var3d(:,:,:,:) + integer ix, k + + nt=nt+1 + do k=lbound(var_block,2),ubound(var_block,2) + do ix=1,size(var_block,1) + var3d(ii1(ix),jj1(ix),k,nt) = var_block(ix,k) + enddo + enddo + end subroutine copy_from_GFS_Data_3d_phys2phys + + pure subroutine copy_from_GFS_Data_2d_int2phys(ii1,jj1,isc,jsc,nt,var2d,var_block) + implicit none + integer, intent(in) :: ii1(:), jj1(:), isc, jsc, var_block(:) + integer, intent(inout) :: nt + real(kind=kind_phys), intent(out) :: var2d(:,:,:) + integer ix + + nt=nt+1 + do ix=1,size(var_block) + var2d(ii1(ix),jj1(ix),nt) = var_block(ix) + enddo + end subroutine copy_from_GFS_Data_2d_int2phys + + pure subroutine copy_from_GFS_Data_2d_stack_phys2phys(ii1,jj1,isc,jsc,nt,var3d,var_block) + ! For copying phy_f2d and phy_fctd + implicit none + integer, intent(in) :: ii1(:), jj1(:), isc, jsc + integer, intent(inout) :: nt + real(kind=kind_phys), intent(in) :: var_block(:,:) + real(kind=kind_phys), intent(out) :: var3d(:,:,:) + integer ix, k + + nt=nt+1 + do k=lbound(var_block,2),ubound(var_block,2) + do ix=1,size(var_block,1) + var3d(ii1(ix),jj1(ix),nt) = var_block(ix,k) + enddo + enddo + end subroutine copy_from_GFS_Data_2d_stack_phys2phys + + pure subroutine copy_from_GFS_Data_3d_int2phys(ii1,jj1,isc,jsc,nt,var3d,var_block) + implicit none + integer, intent(in) :: ii1(:), jj1(:), var_block(:,:), isc, jsc + integer, intent(inout) :: nt + real(kind=kind_phys), intent(out) :: var3d(:,:,:,:) + integer ix, k + + nt=nt+1 + do k=lbound(var_block,2),ubound(var_block,2) + do ix=1,size(var_block,1) + var3d(ii1(ix),jj1(ix),k,nt) = real(var_block(ix,k),kind_phys) + enddo + enddo + end subroutine copy_from_GFS_Data_3d_int2phys + +end module FV3GFS_restart_io_mod diff --git a/io/module_write_internal_state.F90 b/io/module_write_internal_state.F90 index 28e41b1cf..ef7e69f6b 100644 --- a/io/module_write_internal_state.F90 +++ b/io/module_write_internal_state.F90 @@ -48,10 +48,7 @@ module write_internal_state !*** file bundle for output !-------------------------- integer :: FBCount - integer,dimension(:), allocatable :: ncount_attribs - integer,dimension(:), allocatable :: ncount_fields character(128),dimension(:),allocatable :: wrtFB_names - character(128),dimension(:,:),allocatable :: field_names ! !----------------------------------------------------------------------- !*** THE OUTPUT FILE diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 7bf85686e..380ea5975 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -122,7 +122,7 @@ subroutine write_netcdf(wrtfb, filename, & call ESMF_FieldGet(fcstField(i), dimCount=fieldDimCount, array=array, rc=rc); ESMF_ERR_RETURN(rc) if (fieldDimCount > 3) then - write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" + if (mype==0) write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" call ESMF_Finalize(endflag=ESMF_END_ABORT) end if @@ -200,11 +200,11 @@ subroutine write_netcdf(wrtfb, filename, & if (par) then ncerr = nf90_create(trim(filename),& - cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),& + cmode=IOR(NF90_CLOBBER,NF90_NETCDF4),& comm=mpi_comm, info = MPI_INFO_NULL, ncid=ncid); NC_ERR_STOP(ncerr) else ncerr = nf90_create(trim(filename),& - cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),& + cmode=IOR(NF90_CLOBBER,NF90_NETCDF4),& ncid=ncid); NC_ERR_STOP(ncerr) end if @@ -216,13 +216,13 @@ subroutine write_netcdf(wrtfb, filename, & ncerr = nf90_def_dim(ncid, "grid_yt", jm, jm_dimid); NC_ERR_STOP(ncerr) ncerr = nf90_def_dim(ncid, "nchars", 20, ch_dimid); NC_ERR_STOP(ncerr) if (lm > 1) then - call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, rc) - call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, rc) + call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, mype, rc) + call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, mype, rc) end if if (is_cubed_sphere) then ncerr = nf90_def_dim(ncid, "tile", tileCount, tile_dimid); NC_ERR_STOP(ncerr) end if - call add_dim(ncid, "time", time_dimid, wrtgrid, rc) + call add_dim(ncid, "time", time_dimid, wrtgrid, mype, rc) ! define coordinate variables ncerr = nf90_def_var(ncid, "grid_xt", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr) @@ -291,7 +291,7 @@ subroutine write_netcdf(wrtfb, filename, & end if - call get_global_attr(wrtfb, ncid, rc) + call get_global_attr(wrtfb, ncid, mype, rc) ! define variables (fields) @@ -344,7 +344,7 @@ subroutine write_netcdf(wrtfb, filename, & ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & dimids_2d, varids(i)); NC_ERR_STOP(ncerr) else - write(0,*)'Unsupported typekind ', typekind + if (mype==0) write(0,*)'Unsupported typekind ', typekind call ESMF_Finalize(endflag=ESMF_END_ABORT) end if else if (rank == 3) then @@ -384,11 +384,11 @@ subroutine write_netcdf(wrtfb, filename, & ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & dimids_3d, varids(i)); NC_ERR_STOP(ncerr) else - write(0,*)'Unsupported typekind ', typekind + if (mype==0) write(0,*)'Unsupported typekind ', typekind call ESMF_Finalize(endflag=ESMF_END_ABORT) end if else - write(0,*)'Unsupported rank ', rank + if (mype==0) write(0,*)'Unsupported rank ', rank call ESMF_Finalize(endflag=ESMF_END_ABORT) end if if (par) then @@ -513,7 +513,7 @@ subroutine write_netcdf(wrtfb, filename, & end do ncerr = nf90_put_var(ncid, im_varid, values=x); NC_ERR_STOP(ncerr) else - write(0,*)'unknown output_grid ', trim(output_grid(grid_id)) + if (mype==0) write(0,*)'unknown output_grid ', trim(output_grid(grid_id)) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if end if @@ -564,7 +564,7 @@ subroutine write_netcdf(wrtfb, filename, & end do ncerr = nf90_put_var(ncid, jm_varid, values=y); NC_ERR_STOP(ncerr) else - write(0,*)'unknown output_grid ', trim(output_grid(grid_id)) + if (mype==0) write(0,*)'unknown output_grid ', trim(output_grid(grid_id)) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if end if @@ -708,7 +708,7 @@ subroutine write_netcdf(wrtfb, filename, & else - write(0,*)'Unsupported rank ', rank + if (mype==0) write(0,*)'Unsupported rank ', rank call ESMF_Finalize(endflag=ESMF_END_ABORT) end if ! end rank @@ -755,9 +755,10 @@ subroutine write_netcdf(wrtfb, filename, & end subroutine write_netcdf !---------------------------------------------------------------------------------------- - subroutine get_global_attr(fldbundle, ncid, rc) + subroutine get_global_attr(fldbundle, ncid, mype, rc) type(ESMF_FieldBundle), intent(in) :: fldbundle integer, intent(in) :: ncid + integer, intent(in) :: mype integer, intent(out) :: rc ! local variable @@ -766,7 +767,8 @@ subroutine get_global_attr(fldbundle, ncid, rc) character(len=ESMF_MAXSTR) :: attName type(ESMF_TypeKind_Flag) :: typekind - integer :: varival + integer(ESMF_KIND_I4) :: varival_i4 + integer(ESMF_KIND_I8) :: varival_i8 real(ESMF_KIND_R4), dimension(:), allocatable :: varr4list real(ESMF_KIND_R8), dimension(:), allocatable :: varr8list integer :: itemCount @@ -784,8 +786,13 @@ subroutine get_global_attr(fldbundle, ncid, rc) if (typekind==ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & - name=trim(attName), value=varival, rc=rc); ESMF_ERR_RETURN(rc) - ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), varival); NC_ERR_STOP(ncerr) + name=trim(attname), value=varival_i4, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_att(ncid, nf90_global, trim(attname), varival_i4); NC_ERR_STOP(ncerr) + + else if (typekind==ESMF_TYPEKIND_I8) then + call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & + name=trim(attname), value=varival_i8, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_att(ncid, nf90_global, trim(attname), varival_i8); NC_ERR_STOP(ncerr) else if (typekind==ESMF_TYPEKIND_R4) then allocate (varr4list(itemCount)) @@ -806,6 +813,10 @@ subroutine get_global_attr(fldbundle, ncid, rc) name=trim(attName), value=varcval, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), trim(varcval)); NC_ERR_STOP(ncerr) + else + + if (mype==0) write(0,*)'Unsupported typekind ', typekind + call ESMF_Finalize(endflag=ESMF_END_ABORT) end if end do @@ -877,11 +888,12 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) end subroutine get_grid_attr !---------------------------------------------------------------------------------------- - subroutine add_dim(ncid, dim_name, dimid, grid, rc) + subroutine add_dim(ncid, dim_name, dimid, grid, mype, rc) integer, intent(in) :: ncid character(len=*), intent(in) :: dim_name integer, intent(inout) :: dimid type(ESMF_Grid), intent(in) :: grid + integer, intent(in) :: mype integer, intent(out) :: rc ! local variable @@ -927,7 +939,7 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR4) else - write(0,*)'Error in module_write_netcdf.F90(add_dim) unknown typekind for ',trim(dim_name) + if (mype==0) write(0,*)'Error in module_write_netcdf.F90(add_dim) unknown typekind for ',trim(dim_name) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if if (par) then diff --git a/io/module_write_restart_netcdf.F90 b/io/module_write_restart_netcdf.F90 new file mode 100644 index 000000000..259079bb2 --- /dev/null +++ b/io/module_write_restart_netcdf.F90 @@ -0,0 +1,587 @@ +#define ESMF_ERR_RETURN(rc) \ + if (ESMF_LogFoundError(rc, msg="Breaking out of subroutine", line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) + +#define NC_ERR_STOP(status) \ + if (status /= nf90_noerr) write(0,*) "file: ", __FILE__, " line: ", __LINE__, trim(nf90_strerror(status)); \ + if (status /= nf90_noerr) call ESMF_Finalize(endflag=ESMF_END_ABORT) + +module module_write_restart_netcdf + + use esmf + use fms + use netcdf + use mpi + + implicit none + private + public write_restart_netcdf + + logical :: par + + contains + +!---------------------------------------------------------------------------------------- + subroutine write_restart_netcdf(wrtfb, filename, & + use_parallel_netcdf, mpi_comm, mype, & + rc) +! + type(ESMF_FieldBundle), intent(in) :: wrtfb + character(*), intent(in) :: filename + logical, intent(in) :: use_parallel_netcdf + integer, intent(in) :: mpi_comm + integer, intent(in) :: mype + integer, optional,intent(out) :: rc +! +!** local vars + integer, parameter :: i8_kind=8 + integer :: i,j,k,t, istart,iend,jstart,jend + integer :: im, jm, lm + + integer, dimension(:), allocatable :: fldlev + + real(ESMF_KIND_R4), dimension(:,:), pointer :: array_r4 + real(ESMF_KIND_R4), dimension(:,:,:), pointer :: array_r4_cube + real(ESMF_KIND_R4), dimension(:,:,:), pointer :: array_r4_3d + real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: array_r4_3d_cube + + real(ESMF_KIND_R8), dimension(:,:), pointer :: array_r8 + real(ESMF_KIND_R8), dimension(:,:,:), pointer :: array_r8_cube + real(ESMF_KIND_R8), dimension(:,:,:), pointer :: array_r8_3d + real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: array_r8_3d_cube + + real(8), dimension(:), allocatable :: x,y + integer :: fieldCount, fieldDimCount, gridDimCount + integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound + integer, dimension(:), allocatable :: start_idx + + type(ESMF_Field), allocatable :: fcstField(:) + type(ESMF_TypeKind_Flag) :: typekind + type(ESMF_StaggerLoc) :: staggerloc + type(ESMF_TypeKind_Flag) :: attTypeKind + type(ESMF_Grid) :: wrtgrid + type(ESMF_Array) :: array + type(ESMF_DistGrid) :: distgrid + + integer :: attCount + character(len=ESMF_MAXSTR) :: attName, fldName + + integer :: varival + real(4) :: varr4val, dataMin, dataMax + real(4), allocatable, dimension(:) :: compress_err + real(8) :: varr8val + character(len=ESMF_MAXSTR) :: varcval + + integer :: ncerr,ierr + integer :: ncid + integer :: oldMode + integer :: dimid + integer :: im_dimid, im_p1_dimid, jm_dimid, jm_p1_dimid, time_dimid + integer :: im_varid, im_p1_varid, jm_varid, jm_p1_varid, time_varid + integer, dimension(:), allocatable :: dimids_2d, dimids_3d + integer, dimension(:), allocatable :: varids, zaxis_dimids + logical shuffle + + logical :: get_im_jm, found_im_jm + integer :: rank, deCount, localDeCount, dimCount, tileCount + integer :: my_tile, start_i, start_j + integer, dimension(:,:), allocatable :: minIndexPDe, maxIndexPDe + integer, dimension(:,:), allocatable :: minIndexPTile, maxIndexPTile + integer, dimension(:), allocatable :: deToTileMap, localDeToDeMap + logical :: do_io + integer :: par_access + + real(ESMF_KIND_R8), allocatable :: valueListr8(:) + logical :: isPresent, thereAreVerticals, is_restart_core, dynamics_restart_file + integer :: udimCount + character(80), allocatable :: udimList(:) + ! character(32), allocatable :: field_checksums(:) + character(32) :: axis_attr_name + character(32) :: field_checksum +! + is_restart_core = .false. + if ( index(trim(filename),'fv_core.res') > 0 ) is_restart_core = .true. + tileCount = 0 + my_tile = 0 + start_i = -10000000 + start_j = -10000000 + + par = use_parallel_netcdf + do_io = par .or. (mype==0) + + call ESMF_FieldBundleGet(wrtfb, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) + + allocate(compress_err(fieldCount)); compress_err=-999. + allocate(fldlev(fieldCount)) ; fldlev = 0 + allocate(fcstField(fieldCount)) + allocate(varids(fieldCount)) + allocate(zaxis_dimids(fieldCount)) + ! allocate(field_checksums(fieldCount)) + + call ESMF_FieldBundleGet(wrtfb, fieldList=fcstField, grid=wrtGrid, & + itemorderflag=ESMF_ITEMORDER_ADDORDER, & + rc=rc); ESMF_ERR_RETURN(rc) + + call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc); ESMF_ERR_RETURN(rc) + + found_im_jm = .false. + do i=1,fieldCount + call ESMF_FieldGet(fcstField(i), dimCount=fieldDimCount, array=array, rc=rc); ESMF_ERR_RETURN(rc) + call ESMF_FieldGet(fcstField(i), name=fldName, rank=rank, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + + if (fieldDimCount > 3) then + if (mype==0) write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! use first field to determine tile number, grid size, start index etc. + if (is_restart_core) then + get_im_jm = trim(fldName) /= 'u' .and. trim(fldName) /= 'v' ! skip staggered fields + else + get_im_jm = (i==1) + end if + if ( .not.found_im_jm .and. get_im_jm) then + call ESMF_ArrayGet(array, & + distgrid=distgrid, & + dimCount=dimCount, & + deCount=deCount, & + localDeCount=localDeCount, & + tileCount=tileCount, & + rc=rc); ESMF_ERR_RETURN(rc) + + if (tileCount /= 1) then + if (mype==0) write(0,*)"write_restart_netcdf: Only tileCount==1 fields are supported!" + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + + allocate(minIndexPDe(dimCount,deCount)) + allocate(maxIndexPDe(dimCount,deCount)) + allocate(minIndexPTile(dimCount, tileCount)) + allocate(maxIndexPTile(dimCount, tileCount)) + call ESMF_DistGridGet(distgrid, & + minIndexPDe=minIndexPDe, maxIndexPDe=maxIndexPDe, & + minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, & + rc=rc); ESMF_ERR_RETURN(rc) + + allocate(deToTileMap(deCount)) + allocate(localDeToDeMap(localDeCount)) + call ESMF_ArrayGet(array, & + deToTileMap=deToTileMap, & + localDeToDeMap=localDeToDeMap, & + rc=rc); ESMF_ERR_RETURN(rc) + + my_tile = deToTileMap(localDeToDeMap(1)+1) + im = maxIndexPTile(1,1) - minIndexPTile(1,1) + 1 + jm = maxIndexPTile(2,1) - minIndexPTile(2,1) + 1 + start_i = minIndexPDe(1,localDeToDeMap(1)+1) + start_j = minIndexPDe(2,localDeToDeMap(1)+1) + if (.not. par) then + start_i = 1 + start_j = 1 + end if + found_im_jm = .true. + start_i = mod(start_i, im) ! for cubed sphere multi-tile grid + start_j = mod(start_j, jm) ! for cubed sphere multi-tile grid + + deallocate(minIndexPDe) + deallocate(maxIndexPDe) + deallocate(minIndexPTile) + deallocate(maxIndexPTile) + deallocate(deToTileMap) + deallocate(localDeToDeMap) + end if + + if (fieldDimCount > gridDimCount) then + allocate(ungriddedLBound(fieldDimCount-gridDimCount)) + allocate(ungriddedUBound(fieldDimCount-gridDimCount)) + call ESMF_FieldGet(fcstField(i), & + ungriddedLBound=ungriddedLBound, & + ungriddedUBound=ungriddedUBound, rc=rc); ESMF_ERR_RETURN(rc) + fldlev(i) = ungriddedUBound(fieldDimCount-gridDimCount) - & + ungriddedLBound(fieldDimCount-gridDimCount) + 1 + deallocate(ungriddedLBound) + deallocate(ungriddedUBound) + else if (fieldDimCount == 2) then + fldlev(i) = 1 + end if + + end do + + ! create netcdf file and enter define mode + if (do_io) then + + if (par) then + ncerr = nf90_create(trim(filename),& + cmode=IOR(NF90_CLOBBER,NF90_NETCDF4),& + comm=mpi_comm, info = MPI_INFO_NULL, ncid=ncid); NC_ERR_STOP(ncerr) + else + ncerr = nf90_create(trim(filename),& + ! cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),& + cmode=IOR(NF90_CLOBBER,NF90_NETCDF4),& + ncid=ncid); NC_ERR_STOP(ncerr) + end if + + ! disable auto filling. + ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) + + dynamics_restart_file = index(trim(filename),"fv_") > 0 + + ! Unnecessary naming inconsistency + if (dynamics_restart_file) then + axis_attr_name = "axis" + else + axis_attr_name = "cartesian_axis" + end if + + ! define dimensions [xaxis_1, yaxis_1 ,(zaxis_1,...), Time] + if ( .not.is_restart_core ) then + + ncerr = nf90_def_dim(ncid, "xaxis_1", im, im_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "xaxis_1", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr) + + ncerr = nf90_def_dim(ncid, "yaxis_1", jm, jm_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "yaxis_1", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr) + + else + + ncerr = nf90_def_dim(ncid, "xaxis_1", im, im_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "xaxis_1", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr) + + ncerr = nf90_def_dim(ncid, "xaxis_2", im+1, im_p1_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "xaxis_2", NF90_DOUBLE, im_p1_dimid, im_p1_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_p1_varid, trim(axis_attr_name), "X"); NC_ERR_STOP(ncerr) + + ncerr = nf90_def_dim(ncid, "yaxis_1", jm+1, jm_p1_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "yaxis_1", NF90_DOUBLE, jm_p1_dimid, jm_p1_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_p1_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr) + + ncerr = nf90_def_dim(ncid, "yaxis_2", jm, jm_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "yaxis_2", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, trim(axis_attr_name), "Y"); NC_ERR_STOP(ncerr) + + end if + + ! define ungridded (vertical) coordinate variables + do i=1,fieldCount + call ESMF_FieldGet(fcstField(i), name=fldName, rank=rank, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + call ESMF_AttributeGetAttPack(fcstField(i), convention="NetCDF", purpose="FV3", isPresent=isPresent, rc=rc); ESMF_ERR_RETURN(rc) + if (isPresent) then + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + name="ESMF:ungridded_dim_labels", isPresent=isPresent, & + itemCount=udimCount, rc=rc); ESMF_ERR_RETURN(rc) + if (udimCount>1) then + if (mype==0) write(0,*)'udimCount>1 for ', trim(fldName) + ESMF_ERR_RETURN(-1) + end if + if (udimCount>0 .and. isPresent) then + thereAreVerticals = .true. + allocate(udimList(udimCount)) + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", name="ESMF:ungridded_dim_labels", valueList=udimList, rc=rc); ESMF_ERR_RETURN(rc) + ! loop over all ungridded dimension labels + do k=1, udimCount + call write_out_ungridded_dim_atts_from_field(fcstField(i), trim(udimList(k)), dimid, rc=rc); ESMF_ERR_RETURN(rc) + zaxis_dimids(i) = dimid + enddo + deallocate(udimList) + end if + end if + end do + + ncerr = nf90_def_dim(ncid, "Time", NF90_UNLIMITED, time_dimid); NC_ERR_STOP(ncerr) + ! ncerr = nf90_def_dim(ncid, "Time", 1, time_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "Time", NF90_DOUBLE, time_dimid, time_varid); NC_ERR_STOP(ncerr) + if (par) then + ncerr = nf90_var_par_access(ncid, time_varid, NF90_COLLECTIVE); NC_ERR_STOP(ncerr) + end if + + ! Again, unnecessary naming inconsistency + if (dynamics_restart_file) then + ncerr = nf90_put_att(ncid, time_varid, "cartesian_axis", "T"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, time_varid, "units", "time level"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, time_varid, "long_name", "Time"); NC_ERR_STOP(ncerr) + else + ncerr = nf90_put_att(ncid, time_varid, trim(axis_attr_name), "T"); NC_ERR_STOP(ncerr) + end if + + ncerr = nf90_redef(ncid=ncid) + + ! define variables (fields) + allocate(dimids_2d(3)) + allocate(dimids_3d(4)) + + do i=1, fieldCount + call ESMF_FieldGet(fcstField(i), name=fldName, rank=rank, typekind=typekind, staggerloc=staggerloc, rc=rc); ESMF_ERR_RETURN(rc) + + ! FIXME remove this once ESMF_FieldGet above returns correct staggerloc + staggerloc = ESMF_STAGGERLOC_CENTER + if (is_restart_core .and. trim(fldName) == 'u' ) staggerloc = ESMF_STAGGERLOC_EDGE2 + if (is_restart_core .and. trim(fldName) == 'v' ) staggerloc = ESMF_STAGGERLOC_EDGE1 + + ! par_access = NF90_INDEPENDENT + par_access = NF90_COLLECTIVE ! because of time unlimited + ! define variables + if (rank == 2) then + dimids_2d = [im_dimid,jm_dimid, time_dimid] + if (typekind == ESMF_TYPEKIND_R4) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, dimids_2d, varids(i)); NC_ERR_STOP(ncerr) + else if (typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, dimids_2d, varids(i)); NC_ERR_STOP(ncerr) + else + if (mype==0) write(0,*)'Unsupported typekind ', typekind + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + else if (rank == 3) then + if ( .not.is_restart_core ) then + dimids_3d = [im_dimid,jm_dimid,zaxis_dimids(i),time_dimid] + else + if (staggerloc == ESMF_STAGGERLOC_CENTER) then + dimids_3d = [im_dimid,jm_dimid,zaxis_dimids(i),time_dimid] + else if (staggerloc == ESMF_STAGGERLOC_EDGE1) then ! east + dimids_3d = [im_p1_dimid,jm_dimid,zaxis_dimids(i),time_dimid] + else if (staggerloc == ESMF_STAGGERLOC_EDGE2) then ! south + dimids_3d = [im_dimid,jm_p1_dimid,zaxis_dimids(i),time_dimid] + else + if (mype==0) write(0,*)'Unsupported staggerloc ', staggerloc + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + if (typekind == ESMF_TYPEKIND_R4) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, dimids_3d, varids(i)); NC_ERR_STOP(ncerr) + else if (typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, dimids_3d, varids(i)); NC_ERR_STOP(ncerr) + else + if (mype==0) write(0,*)'Unsupported typekind ', typekind + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + else + if (mype==0) write(0,*)'Unsupported rank ', rank + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + if (par) then + ncerr = nf90_var_par_access(ncid, varids(i), par_access); NC_ERR_STOP(ncerr) + end if + + end do ! i=1,fieldCount + + ncerr = nf90_put_att(ncid, NF90_GLOBAL, "NumFilesInSet", 1); NC_ERR_STOP(ncerr) + +! end of define mode + ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr) + + if (allocated(valueListr8)) deallocate(valueListr8) + allocate (valueListr8(im)) + valueListr8 = (/(i, i=1,im)/) + ncerr = nf90_put_var(ncid, im_varid, values=valueListr8); NC_ERR_STOP(ncerr) + + if (allocated(valueListr8)) deallocate(valueListr8) + allocate (valueListr8(jm)) + valueListr8 = (/(i, i=1,jm)/) + ncerr = nf90_put_var(ncid, jm_varid, values=valueListr8); NC_ERR_STOP(ncerr) + + if ( is_restart_core ) then + if (allocated(valueListr8)) deallocate(valueListr8) + allocate (valueListr8(im+1)) + valueListr8 = (/(i, i=1,im+1)/) + ncerr = nf90_put_var(ncid, im_p1_varid, values=valueListr8); NC_ERR_STOP(ncerr) + + if (allocated(valueListr8)) deallocate(valueListr8) + allocate (valueListr8(jm+1)) + valueListr8 = (/(i, i=1,jm+1)/) + ncerr = nf90_put_var(ncid, jm_p1_varid, values=valueListr8); NC_ERR_STOP(ncerr) + end if + + ncerr = nf90_put_var(ncid, time_varid, values=[1]); NC_ERR_STOP(ncerr) + end if + + ! write variables (fields) + do i=1, fieldCount + + call ESMF_FieldGet(fcstField(i),name=fldName,rank=rank,typekind=typekind,staggerloc=staggerloc, rc=rc); ESMF_ERR_RETURN(rc) + + ! FIXME remove this once ESMF_FieldGet above returns correct staggerloc + staggerloc = ESMF_STAGGERLOC_CENTER + if (is_restart_core .and. trim(fldName) == 'u' ) staggerloc = ESMF_STAGGERLOC_EDGE2 + if (is_restart_core .and. trim(fldName) == 'v' ) staggerloc = ESMF_STAGGERLOC_EDGE1 + + if (rank == 2) then + + if (allocated(start_idx)) deallocate(start_idx) + allocate(start_idx(3)) + start_idx = [start_i,start_j, 1] + + if (typekind == ESMF_TYPEKIND_R4) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r4, rc=rc); ESMF_ERR_RETURN(rc) + write(field_checksum,'(Z16)') mpp_chksum(array_r4) + if (par) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4, start=start_idx); NC_ERR_STOP(ncerr) + else + allocate(array_r4(im,jm)) + call ESMF_FieldGather(fcstField(i), array_r4, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4, start=start_idx); NC_ERR_STOP(ncerr) + end if + deallocate(array_r4) + end if + else if (typekind == ESMF_TYPEKIND_R8) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r8, rc=rc); ESMF_ERR_RETURN(rc) + write(field_checksum,'(Z16)') mpp_chksum(array_r8) + if (par) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + else + allocate(array_r8(im,jm)) + call ESMF_FieldGather(fcstField(i), array_r8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + end if + deallocate(array_r8) + end if + end if + + else if (rank == 3) then + + if (allocated(start_idx)) deallocate(start_idx) + allocate(start_idx(4)) + start_idx = [start_i,start_j,1, 1] + + if (typekind == ESMF_TYPEKIND_R4) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r4_3d, rc=rc); ESMF_ERR_RETURN(rc) + write(field_checksum,'(Z16)') mpp_chksum(array_r4_3d) + if (par) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d, start=start_idx); NC_ERR_STOP(ncerr) + else + if (staggerloc == ESMF_STAGGERLOC_CENTER) then + allocate(array_r4_3d(im,jm,fldlev(i))) + else if (staggerloc == ESMF_STAGGERLOC_EDGE1) then ! east + allocate(array_r4_3d(im+1,jm,fldlev(i))) + else if (staggerloc == ESMF_STAGGERLOC_EDGE2) then ! south + allocate(array_r4_3d(im,jm+1,fldlev(i))) + else + if (mype==0) write(0,*)'Unsupported staggerloc ', staggerloc + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + call ESMF_FieldGather(fcstField(i), array_r4_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (mype==0) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d, start=start_idx); NC_ERR_STOP(ncerr) + end if + deallocate(array_r4_3d) + end if + else if (typekind == ESMF_TYPEKIND_R8) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r8_3d, rc=rc); ESMF_ERR_RETURN(rc) + write(field_checksum,'(Z16)') mpp_chksum(array_r8_3d) + if (par) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_3d, start=start_idx); NC_ERR_STOP(ncerr) + else + if (staggerloc == ESMF_STAGGERLOC_CENTER) then + allocate(array_r8_3d(im,jm,fldlev(i))) + else if (staggerloc == ESMF_STAGGERLOC_EDGE1) then ! east + allocate(array_r8_3d(im+1,jm,fldlev(i))) + else if (staggerloc == ESMF_STAGGERLOC_EDGE2) then ! south + allocate(array_r8_3d(im,jm+1,fldlev(i))) + else + if (mype==0) write(0,*)'Unsupported staggerloc ', staggerloc + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + call ESMF_FieldGather(fcstField(i), array_r8_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (mype==0) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_3d, start=start_idx); NC_ERR_STOP(ncerr) + end if + deallocate(array_r8_3d) + end if + end if ! end typekind + + else + + if (mype==0) write(0,*)'Unsupported rank ', rank + call ESMF_Finalize(endflag=ESMF_END_ABORT) + + end if ! end rank + if (do_io) then + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'checksum', field_checksum); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + end if + + end do ! end fieldCount + + if (do_io) then + ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) + end if + +contains + + subroutine write_out_ungridded_dim_atts_from_field(field, dimLabel, dimid, rc) + + type(ESMF_Field),intent(in) :: field + character(len=*),intent(in) :: dimLabel + integer, intent(out) :: dimid + integer, intent(out) :: rc + + real(ESMF_KIND_R4), allocatable :: valueListr4(:) + real(ESMF_KIND_R8), allocatable :: valueListr8(:) + integer :: valueCount, udimCount + integer :: ncerr, varid, ind + integer :: itemCount, attCount + character(len=80), allocatable :: attNameList(:) + character(len=80) :: attName + type(ESMF_TypeKind_Flag) :: typekind + character(len=80) :: valueS + integer :: valueI4 + real(ESMF_KIND_R4) :: valueR4 + real(ESMF_KIND_R8) :: valueR8 + + ! inquire if NetCDF file already contains this ungridded dimension variable + ncerr = nf90_inq_varid(ncid, trim(dimLabel), varid=varid) + if (ncerr == NF90_NOERR) then + ! if it does inquire dimid, it must be defined already + ncerr = nf90_inq_dimid(ncid, trim(dimLabel), dimid=dimid) + if (ncerr == NF90_NOERR) then + return + else + if (mype==0) write(0,*) 'in write_out_ungridded_dim_atts: ERROR missing dimid for already defined ungridded dimension variable ' + ESMF_ERR_RETURN(-1) + endif + endif + + ! the variable does not exist in the NetCDF file yet -> add it + ! access the undistributed dimension attribute on the grid + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", name=trim(dimLabel), itemCount=valueCount, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + if( typekind == ESMF_TYPEKIND_R4 ) then + allocate(valueListr4(valueCount)) + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", name=trim(dimLabel), valueList=valueListr4, rc=rc); ESMF_ERR_RETURN(rc) + else if ( typekind == ESMF_TYPEKIND_R8) then + allocate(valueListr8(valueCount)) + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", name=trim(dimLabel), valueList=valueListr8, rc=rc); ESMF_ERR_RETURN(rc) + else + if (mype==0) write(0,*) 'in write_out_ungridded_dim_atts: ERROR unknown typekind' + ESMF_ERR_RETURN(-1) + endif + + ! now add it to the NetCDF file + ncerr = nf90_inq_dimid(ncid, trim(dimLabel), dimid=dimid) + if (ncerr /= NF90_NOERR) then + ! dimension does not yet exist, and must be defined + ncerr = nf90_def_dim(ncid, trim(dimLabel), valueCount, dimid=dimid); NC_ERR_STOP(ncerr); NC_ERR_STOP(ncerr) + endif + if( typekind == ESMF_TYPEKIND_R4 ) then + ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_FLOAT, dimids=(/dimid/), varid=varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varid, trim(axis_attr_name), "Z"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, varid, values=valueListr4); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + deallocate(valueListr4) + else if(typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_DOUBLE, dimids=(/dimid/), varid=varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varid, trim(axis_attr_name), "Z"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, varid, values=valueListr8); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + deallocate(valueListr8) + endif + end subroutine write_out_ungridded_dim_atts_from_field + + end subroutine write_restart_netcdf + +!---------------------------------------------------------------------------------------- +end module module_write_restart_netcdf diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index c3b76701a..ba4d44e16 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -28,6 +28,7 @@ module module_wrt_grid_comp ! use mpi use esmf + use fms use write_internal_state use module_fv3_io_def, only : num_pes_fcst, & @@ -41,6 +42,7 @@ module module_wrt_grid_comp stdlat1, stdlat2, dx, dy, iau_offset, & ideflate, lflname_fulltime use module_write_netcdf, only : write_netcdf + use module_write_restart_netcdf use physcons, only : pi => con_pi #ifdef INLINE_POST use post_fv3, only : post_run_fv3 @@ -67,11 +69,13 @@ module module_wrt_grid_comp integer,save :: idate(7) logical,save :: write_nsflip logical,save :: change_wrtidate=.false. + integer,save :: frestart(999) = -1 + logical :: lprnt ! !----------------------------------------------------------------------- ! type(ESMF_FieldBundle) :: gridFB - integer :: FBcount + integer :: FBCount character(len=esmf_maxstr),allocatable :: fcstItemNameList(:) ! !----------------------------------------------------------------------- @@ -158,6 +162,7 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, type(ESMF_DELayout) :: delayout type(ESMF_Grid) :: fcstGrid type(ESMF_Grid), allocatable :: wrtGrid(:) + type(ESMF_Grid) :: actualWrtGrid type(ESMF_Array) :: array type(ESMF_Field) :: field_work, field type(ESMF_Decomp_Flag) :: decompflagPTile(2,6) @@ -172,6 +177,7 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, integer, allocatable :: gridToFieldMap(:) integer, allocatable :: ungriddedLBound(:) integer, allocatable :: ungriddedUBound(:) + type(ESMF_StaggerLoc) :: staggerloc character(len=80) :: attName character(len=80), allocatable :: attNameList(:),attNameList2(:) type(ESMF_TypeKind_Flag), allocatable :: typekindList(:) @@ -180,6 +186,7 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, real(ESMF_KIND_R4) :: valueR4 real(ESMF_KIND_R8) :: valueR8 logical, allocatable :: is_moving(:) + logical :: isPresent integer :: attCount, jidx, idx, noutfile character(19) :: newdate @@ -197,11 +204,10 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, character(256) :: cf_open, cf_close character(256) :: gridfile integer :: num_output_file -! - logical :: lprnt - integer :: grid_id - logical :: top_parent_is_global + type(ESMF_DistGrid) :: acceptorDG, newAcceptorDG + integer :: grid_id + logical :: top_parent_is_global ! !----------------------------------------------------------------------- !*********************************************************************** @@ -239,6 +245,9 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, last_write_task = ntasks -1 lprnt = lead_write_task == wrt_int_state%mype + call fms_init(wrt_mpi_comm) + call mpp_init() + ! print *,'in wrt, lead_write_task=', & ! lead_write_task,'last_write_task=',last_write_task, & ! 'mype=',wrt_int_state%mype,'jidx=',jidx,' comm=',wrt_mpi_comm @@ -777,18 +786,18 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, ,startTime=wrt_int_state%IO_BASETIME & !<-- The Clock's starting time ,rc =RC) - call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), & - m=idate(5),s=idate(6),rc=rc) + call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3), & + h=idate(4), m=idate(5), s=idate(6),rc=rc) ! if (lprnt) write(0,*) 'in wrt initial, io_baseline time=',idate,'rc=',rc idate(7) = 1 wrt_int_state%idate = idate wrt_int_state%fdate = idate ! update IO-BASETIME and idate on write grid comp when IAU is enabled - if(iau_offset > 0 ) then + if (iau_offset > 0) then call ESMF_TimeIntervalSet(IAU_offsetTI, h=iau_offset, rc=rc) wrt_int_state%IO_BASETIME = wrt_int_state%IO_BASETIME + IAU_offsetTI - call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3),h=idate(4), & - m=idate(5),s=idate(6),rc=rc) + call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=idate(1),mm=idate(2),dd=idate(3), & + h=idate(4), m=idate(5), s=idate(6),rc=rc) wrt_int_state%idate = idate change_wrtidate = .true. if (lprnt) print *,'in wrt initial, with iau, io_baseline time=',idate,'rc=',rc @@ -797,19 +806,13 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, !--- Look at the incoming FieldBundles in the imp_state_write, and mirror them as 'output_' bundles ! call ESMF_StateGet(imp_state_write, itemCount=FBCount, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - wrt_int_state%FBCount = FBCount - -! if (lprnt) write(0,*) 'in wrt,fcst FBCount=',FBCount -! if (lprnt) write(0,*) 'in wrt,fcst wrt_int_state%FBCount=',wrt_int_state%FBCount + ! if (lprnt) write(0,*)'wrt_initialize_p1: FBCount=',FBCount, ' from imp_state_write' allocate(fcstItemNameList(FBCount), fcstItemTypeList(FBCount)) - allocate(wrt_int_state%wrtFB_names(FBCount)) - allocate(wrt_int_state%ncount_fields(FBCount),wrt_int_state%ncount_attribs(FBCount)) - allocate(wrt_int_state%field_names(2000,FBCount)) - allocate(outfilename(2000,FBcount)) + allocate(wrt_int_state%wrtFB_names(FBCount)) ! this array should be allocated as wrt_int_state%FBCount long not FBCount + allocate(outfilename(2000,FBCount)) outfilename = '' call ESMF_StateGet(imp_state_write, itemNameList=fcstItemNameList, & @@ -826,16 +829,25 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, call ESMF_StateGet(imp_state_write, itemName=fcstItemNameList(i), & fieldbundle=fcstFB, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(fcstFB, convention="NetCDF", purpose="FV3", & - name="grid_id", value=grid_id, rc=rc) + name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return -!--- get grid dim count - call ESMF_GridGet(wrtGrid(grid_id), dimCount=gridDimCount, rc=rc) + call ESMF_AttributeGet(fcstFB, convention="NetCDF", purpose="FV3-nooutput", & + name="frestart", valueList=frestart, isPresent=isPresent, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (isPresent) then + ! if (lprnt) write(0,*)'wrt_initialize_p1: frestart(1:10) = ',frestart(1:10) + call ESMF_AttributeRemove(fcstFB, convention="NetCDF", purpose="FV3-nooutput", name="frestart", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + endif + + +!--- get grid dim count + ! call ESMF_GridGet(wrtGrid(grid_id), dimCount=gridDimCount, rc=rc) + ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! create a mirrored 'output_' FieldBundle and add it to importState fieldbundle = ESMF_FieldBundleCreate(name="output_"//trim(fcstItemNameList(i)), rc=rc) @@ -865,24 +877,39 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, endif ! deal with all of the Fields inside this fcstFB - call ESMF_FieldBundleGet(fcstFB, fieldCount=fieldCount, rc=rc) + call ESMF_FieldBundleGet(fcstFB, fieldCount=fieldCount, grid=fcstGrid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (fieldCount > 0) then - wrt_int_state%ncount_fields(i) = fieldCount - allocate(fcstField(fieldCount)) call ESMF_FieldBundleGet(fcstFB, fieldList=fcstField, & itemorderflag=ESMF_ITEMORDER_ADDORDER, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + actualWrtGrid = wrtGrid(grid_id) + + ! If this is a 'restart' bundle the actual grid that the output field ('field_work' below) is created on + ! must be the same grid as forecast grid, not the output grid for this grid_id (wrtGrid(grid_id)). + ! For 'cubed_sphere_grid' these are the same, but for all other output grids (like Lambert) they are not. + if (fcstItemNameList(i)(1:8) == 'restart_') then + ! create a grid from fcstGrid on forecast grid comp, by rebalancing distgrid to the local PETs + ! access the acceptor DistGrid + call ESMF_GridGet(fcstGrid, distgrid=acceptorDG, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! rebalance the acceptor DistGrid across the local PETs + newAcceptorDG = ESMF_DistGridCreate(acceptorDG, balanceflag=.true., rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + actualWrtGrid = ESMF_GridCreate(fcstGrid, newAcceptorDG, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + end if ! end of setting actualWrtGrid for restart bundle + do j=1, fieldCount - call ESMF_FieldGet(fcstField(j), typekind=typekind, & - dimCount=fieldDimCount, name=fieldName, grid=fcstGrid, rc=rc) + call ESMF_FieldGet(fcstField(j), typekind=typekind, dimCount=fieldDimCount, name=fieldName, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_GridGet(actualWrtGrid, dimCount=gridDimCount, rc=rc) ! use actualWrtGrid instead of wrtGrid(grid_id) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(gridToFieldMap(gridDimCount)) @@ -891,38 +918,30 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, call ESMF_FieldGet(fcstField(j), gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound, & - rc=rc) - call ESMF_LogWrite("after field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) + staggerloc=staggerloc, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! if (lprnt) print *,'in wrt,fcstfld,fieldname=', & ! trim(fieldname),'fieldDimCount=',fieldDimCount,'gridDimCount=',gridDimCount, & ! 'gridToFieldMap=',gridToFieldMap,'ungriddedLBound=',ungriddedLBound, & ! 'ungriddedUBound=',ungriddedUBound,'rc=',rc - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - -! create the output field - - call ESMF_LogWrite("call field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) - field_work = ESMF_FieldCreate(wrtGrid(grid_id), typekind, name=fieldName, & +! create the output field on output grid + field_work = ESMF_FieldCreate(actualWrtGrid, typekind, name=fieldName, & ! use actualWrtGrid instead of wrtGrid(grid_id) + staggerloc=staggerloc, & gridToFieldMap=gridToFieldMap, & ungriddedLBound=ungriddedLBound, & ungriddedUBound=ungriddedUBound, rc=rc) - call ESMF_LogWrite("aft call field create on wrt comp",ESMF_LOGMSG_INFO,rc=RC) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - wrt_int_state%field_names(j,i) = trim(fieldName) - call ESMF_AttributeCopy(fcstField(j), field_work, attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return -! + ! get output file name call ESMF_AttributeGet(fcstField(j), convention="NetCDF", purpose="FV3", & name="output_file", value=outfile_name, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_LogWrite("bf fcstfield, get output_file "//trim(outfile_name)//" "//trim(fieldName),ESMF_LOGMSG_INFO,rc=RC) if (trim(outfile_name) /= '') then outfilename(j,i) = trim(outfile_name) @@ -964,7 +983,8 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, deallocate(gridToFieldMap, ungriddedLBound, ungriddedUBound) enddo ! - call ESMF_AttributeCopy(fcstGrid, wrtGrid(grid_id), & + ! call ESMF_AttributeCopy(fcstGrid, wrtGrid(grid_id), & + call ESMF_AttributeCopy(fcstGrid, actualWrtGrid , & attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -978,42 +998,52 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, return endif - !end FBCount enddo ! !loop over all items in the imp_state_write and count output FieldBundles - call get_outfile(FBcount, outfilename,FBlist_outfilename,noutfile) + ! if (lprnt) write(0,*)'wrt_initialize_p1: before get_outfile FBCount =', FBCount + call get_outfile(FBCount, outfilename, FBlist_outfilename, noutfile) wrt_int_state%FBCount = noutfile + ! if (lprnt) write(0,*)'wrt_initialize_p1: wrt_int_state%FBCount = noutfile ', wrt_int_state%FBCount, noutfile ! !create output field bundles - allocate(wrt_int_state%wrtFB(wrt_int_state%FBcount)) - do i=1, wrt_int_state%FBcount + allocate(wrt_int_state%wrtFB(wrt_int_state%FBCount)) + ! if (lprnt) write(0,*)'wrt_initialize_p1: allocated ',wrt_int_state%FBCount, ' wrt_int_state%wrtFB' + + do i=1, wrt_int_state%FBCount wrt_int_state%wrtFB_names(i) = trim(FBlist_outfilename(i)) wrt_int_state%wrtFB(i) = ESMF_FieldBundleCreate(name=trim(wrt_int_state%wrtFB_names(i)), rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! if (lprnt) write(0,*)'wrt_initialize_p1: created wrtFB ',i, ' with name ', trim(wrt_int_state%wrtFB_names(i)) - do n=1, FBcount + ! if (lprnt) write(0,*)'wrt_initialize_p1: loop over ', FBCount, ' forecast bundles' + do n=1, FBCount call ESMF_StateGet(imp_state_write, itemName="output_"//trim(fcstItemNameList(n)), & fieldbundle=fcstFB, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ! if (lprnt) write(0,*)'wrt_initialize_p1: got forecast bundle ', "output_"//trim(fcstItemNameList(n)) + ! if (lprnt) write(0,*)'wrt_initialize_p1: is ', trim(fcstItemNameList(n)), ' == ', trim(FBlist_outfilename(i)) if( index(trim(fcstItemNameList(n)),trim(FBlist_outfilename(i))) == 1 ) then ! ! copy the fcstfield bundle Attributes to the output field bundle + ! if (lprnt) write(0,*)'wrt_initialize_p1: copy atts/fields from ', "output_"//trim(fcstItemNameList(n)), ' to ', trim(wrt_int_state%wrtFB_names(i)) call ESMF_AttributeCopy(fcstFB, wrt_int_state%wrtFB(i), & attcopy=ESMF_ATTCOPY_REFERENCE, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_AttributeGet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", & - name="grid_id", value=grid_id, rc=rc) + name="grid_id", value=grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return -! + ! if (lprnt) write(0,*)'wrt_initialize_p1: got grid_id for wrtFB ', i, ' grid_id =', grid_id, trim(output_grid(grid_id)) + call ESMF_FieldBundleGet(fcstFB, fieldCount=fieldCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -1044,7 +1074,7 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, enddo deallocate(fcstField, fieldnamelist) - endif + endif ! index(trim(fcstItemNameList(n)),trim(FBlist_outfilename(i))) ! add output grid related attributes @@ -1163,8 +1193,8 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock, end if - enddo ! end FBcount - enddo ! end wrt_int_state%FBcount + enddo ! end FBCount + enddo ! end wrt_int_state%FBCount ! ! add time Attribute ! look at the importState attributes and copy those starting with "time" @@ -1619,7 +1649,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) type(ESMF_TimeInterval) :: io_currtimediff type(ESMF_Grid) :: fbgrid, wrtGrid type(ESMF_State),save :: stateGridFB - type(optimizeT), save :: optimize(4) + type(optimizeT), save :: optimize(40) ! FIXME type(ESMF_GridComp), save, allocatable :: compsGridFB(:) type(ESMF_RouteHandle) :: rh type(ESMF_RegridMethod_Flag) :: regridmethod @@ -1631,17 +1661,19 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) integer :: i,j,n,mype,nolog, grid_id, localPet ! integer :: nf_hours,nf_seconds,nf_minutes + integer :: fcst_seconds real(ESMF_KIND_R8) :: nfhour ! - integer :: nbdl, date(6), ndig, nnnn + integer :: nbdl, cdate(6), ndig, nnnn integer :: step=1 ! logical :: opened logical :: lmask_fields ! - character(esmf_maxstr) :: filename,compname, traceString + character(esmf_maxstr) :: filename,compname,wrtFBName,traceString character(40) :: cfhour, cform character(20) :: time_iso + character(15) :: time_restart ! type(ESMF_Grid) :: grid type(ESMF_Info) :: info @@ -1667,8 +1699,10 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) real(kind=8) :: tbeg real(kind=8) :: wbeg,wend - logical :: use_parallel_netcdf - logical :: lprnt + logical :: use_parallel_netcdf + real, allocatable :: output_fh(:) + logical :: is_restart_bundle + integer :: tileCount ! !----------------------------------------------------------------------- !*********************************************************************** @@ -1698,26 +1732,26 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) call ESMF_VMGetCurrent(VM,rc=RC) mype = wrt_int_state%mype - lprnt = mype == lead_write_task ! print *,'in wrt run, mype=',mype,'lead_write_task=',lead_write_task + + call ESMF_InfoGetFromHost(imp_state_write, info=info, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_InfoGetAlloc(info, key="output_fh", values=output_fh, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! !----------------------------------------------------------------------- !*** get current time and elapsed forecast time - call ESMF_ClockGet(clock=CLOCK, currTime=CURRTIME, rc=rc) + call ESMF_ClockGet(clock=CLOCK, currTime=currTime, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_TimeGet(time=currTime,yy=date(1),mm=date(2),dd=date(3),h=date(4), & - m=date(5),s=date(6),rc=rc) + call ESMF_TimeGet(time=currTime,yy=cdate(1),mm=cdate(2),dd=cdate(3), & + h=cdate(4), m=cdate(5), s=cdate(6),rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return wrt_int_state%fdate(7) = 1 - wrt_int_state%fdate(1:6) = date(1:6) - write(time_iso,'(I4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2,"Z")') date(1:6) - - call ESMF_TimeGet(time=wrt_int_state%IO_BASETIME,yy=date(1),mm=date(2),dd=date(3),h=date(4), & - m=date(5),s=date(6),rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + wrt_int_state%fdate(1:6) = cdate(1:6) + write(time_iso,'(I4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2,"Z")') cdate(1:6) io_currtimediff = currtime - wrt_int_state%IO_BASETIME @@ -1944,8 +1978,13 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) trim(output_grid(n)) == 'rotated_latlon_moving' .or. & trim(output_grid(n)) == 'lambert_conformal') then - !mask fields according to sfc pressure + !mask fields according to sfc pressure, only history bundles do nbdl=1, wrt_int_state%FBCount + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), name=wrtFBName, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__,file=__FILE__)) return + + if (wrtFBName(1:8) == 'restart_') cycle + call mask_fields(wrt_int_state%wrtFB(nbdl),rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo @@ -1978,10 +2017,32 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) ! ** now loop through output field bundle !----------------------------------------------------------------------- - if ( wrt_int_state%output_history ) then + call ESMF_TimeIntervalGet(timeinterval=io_currtimediff, s=fcst_seconds, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ! fcst_seconds is number of seconds in io_currtimediff, which is time interval between currenttime and io_basetime. + ! io_basetime has been adjusted by iau_offset in initialize phase. + ! Since output_fh and frestart and NOT adjusted by iau_offset, in order to compare + ! them with fcst_seconds, we must also adjust fcst_seconds by iau_offset + if (iau_offset > 0) then + fcst_seconds = fcst_seconds + iau_offset*3600 + endif + + if ( (wrt_int_state%output_history .and. ANY(nint(output_fh(:)*3600.0) == fcst_seconds)) .or. ANY(frestart(:) == fcst_seconds) ) then + ! if (lprnt) write(0,*)'wrt_run: loop over wrt_int_state%FBCount ',wrt_int_state%FBCount, ' nfhour ', nfhour, ' cdate ', cdate(1:6) file_loop_all: do nbdl=1, wrt_int_state%FBCount -! + + ! if (lprnt) write(0,*)'wrt_run: nbdl = ',nbdl, ' fb name ',trim(wrt_int_state%wrtFB_names(nbdl)) + + is_restart_bundle = .false. + if (wrt_int_state%wrtFB_names(nbdl)(1:8) == 'restart_') then + is_restart_bundle = .true. + if (.not.(ANY(frestart(:) == fcst_seconds))) cycle + else + if (.not.(wrt_int_state%output_history .and. ANY(nint(output_fh(:)*3600.0) == fcst_seconds))) cycle + endif + ! get grid_id call ESMF_AttributeGet(wrt_int_state%wrtFB(nbdl), convention="NetCDF", purpose="FV3", & name="grid_id", value=grid_id, rc=rc) @@ -2072,7 +2133,32 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) endif endif - filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nc' + if (is_restart_bundle) then + write(time_restart,'(I4,I2.2,I2.2,".",I2.2,I2.2,I2.2)') cdate(1:6) + + ! strip leading 'restart_' from a bundle name and replace it with a directory name 'RESTART/' to create actual file name + filename = 'RESTART/'//trim(time_restart)//'.'//trim(wrt_int_state%wrtFB_names(nbdl)(9:))//'.nc' + + ! I hate this kind of inconsistencies + ! If it's a restart bundle and the output grid is not cubed sphere and the output restart file is + ! from dycore (ie. fv_core, fv_srf_wnd, fv_tracer) append 'tile1' to the end of the file name. + ! As opposed to physics restart files (phy_data, sfc_data) which do not have 'tile1' appended. + ! Why can't we have consistent naming? + + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), grid=grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_GridGet(grid, tileCount=tileCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (tileCount == 1) then ! non cubed sphere restart bundles + if (wrt_int_state%wrtFB_names(nbdl)(9:11) == 'fv_') then ! 'dynamics' restart bundles, append 'tile1' + filename = 'RESTART/'//trim(time_restart)//'.'//trim(wrt_int_state%wrtFB_names(nbdl)(9:))//'.tile1'//'.nc' + endif + endif + + else ! history bundle + filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nc' + endif if(mype == lead_write_task) print *,'in wrt run,filename= ',nbdl,trim(filename) ! @@ -2106,48 +2192,66 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) call ESMF_Finalize(endflag=ESMF_END_ABORT) endif + wbeg = MPI_Wtime() + + if (is_restart_bundle) then ! restart bundle + ! restart bundles are always on forecast grid, either cubed sphere or regional/nest + + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), grid=grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_GridGet(grid, tileCount=tileCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (tileCount == 6) then ! restart bundle is on cubed sphere + call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & + convention="NetCDF", purpose="FV3", & + status=ESMF_FILESTATUS_REPLACE, & + state=stateGridFB, comps=compsGridFB,rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMFproto_FieldBundleWrite(wrt_int_state%wrtFB(nbdl), & + filename=trim(filename), convention="NetCDF", & + purpose="FV3", status=ESMF_FILESTATUS_OLD, & + timeslice=step, state=optimize(nbdl)%state, & + comps=optimize(nbdl)%comps, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + call write_restart_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .true., wrt_mpi_comm, wrt_int_state%mype, & + rc) + endif ! cubed sphere vs. regional/nest write grid + + else ! history bundle if (trim(output_grid(grid_id)) == 'cubed_sphere_grid') then - wbeg = MPI_Wtime() if (trim(output_file(nnnn)) == 'netcdf_parallel') then call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & .true., wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) else - call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & - convention="NetCDF", purpose="FV3", & - status=ESMF_FILESTATUS_REPLACE, & + call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & + convention="NetCDF", purpose="FV3", & + status=ESMF_FILESTATUS_REPLACE, & state=stateGridFB, comps=compsGridFB,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMFproto_FieldBundleWrite(wrt_int_state%wrtFB(nbdl), & - filename=trim(filename), convention="NetCDF", & - purpose="FV3", status=ESMF_FILESTATUS_OLD, & - timeslice=step, state=optimize(nbdl)%state, & - comps=optimize(nbdl)%comps, rc=rc) + call ESMFproto_FieldBundleWrite(wrt_int_state%wrtFB(nbdl), & + filename=trim(filename), convention="NetCDF", & + purpose="FV3", status=ESMF_FILESTATUS_OLD, & + timeslice=step, state=optimize(nbdl)%state, & + comps=optimize(nbdl)%comps, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return end if - wend = MPI_Wtime() - if (lprnt) then - write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & - ,' at Fcst ',NF_HOURS,':',NF_MINUTES - endif else if (trim(output_grid(grid_id)) == 'gaussian_grid' .or. & trim(output_grid(grid_id)) == 'global_latlon') then - wbeg = MPI_Wtime() call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & use_parallel_netcdf, wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) - wend = MPI_Wtime() - if (lprnt) then - write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & - ,' at Fcst ',NF_HOURS,':',NF_MINUTES - endif else if (trim(output_grid(grid_id)) == 'regional_latlon' .or. & trim(output_grid(grid_id)) == 'regional_latlon_moving' .or. & @@ -2157,13 +2261,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) !mask fields according to sfc pressure if( .not. lmask_fields ) then - wbeg = MPI_Wtime() call mask_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - wend = MPI_Wtime() - if (lprnt) then - write(*,'(A,F10.5,A,I4.2,A,I2.2)')' mask_fields time is ',wend-wbeg - endif endif if (nbits(grid_id) /= 0) then @@ -2171,27 +2270,27 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if - wbeg = MPI_Wtime() call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & use_parallel_netcdf, wrt_mpi_comm,wrt_int_state%mype, & grid_id,rc) - wend = MPI_Wtime() - if (lprnt) then - write(*,'(A15,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(output_file(nnnn)),' Write Time is ',wend-wbeg & - ,' at Fcst ',NF_HOURS,':',NF_MINUTES - endif else ! unknown output_grid call ESMF_LogWrite("wrt_run: Unknown output_grid",ESMF_LOGMSG_ERROR,rc=RC) call ESMF_Finalize(endflag=ESMF_END_ABORT) - endif + endif + endif ! restart or history bundle + wend = MPI_Wtime() + if (lprnt) then + write(*,'(A48,A,F10.5,A,I4.2,A,I2.2,1X,A)')trim(filename),' Write Time is ',wend-wbeg & + ,' at Fcst ',NF_HOURS,':',NF_MINUTES + endif - enddo file_loop_all + enddo file_loop_all ! end output history - endif + endif ! if ( wrt_int_state%output_history ) ! !** write out log file ! @@ -2210,7 +2309,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) write_run_tim = MPI_Wtime() - tbeg ! IF (lprnt) THEN - WRITE(*,'(A,F10.5,A,I4.2,A,I2.2)')' total Write Time is ',write_run_tim & + write(*,'(A48,A,F10.5,A,I4.2,A,I2.2,1X,A)')'------- total',' Write Time is ',write_run_tim & ,' at Fcst ',NF_HOURS,':',NF_MINUTES ENDIF ! @@ -2261,6 +2360,8 @@ subroutine wrt_finalize(wrt_comp, imp_state_write, exp_state_write, clock, rc) if (ESMF_LogFoundDeallocError(statusToCheck=stat, & msg="Deallocation of internal state memory failed.", & line=__LINE__, file=__FILE__)) return + + call fms_end ! !----------------------------------------------------------------------- ! @@ -3052,6 +3153,7 @@ subroutine ioCompSS(comp, rc) subroutine ioCompRun(comp, importState, exportState, clock, rc) use netcdf + type(ESMF_GridComp) :: comp type(ESMF_State) :: importState, exportState type(ESMF_Clock) :: clock @@ -3066,7 +3168,7 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) type(ESMF_FileStatus_Flag) :: status character(len=80) :: itemNameList(1) - integer :: localPet, i, j, k, ind + integer :: localPet, petCount, i, j, k, ind type(ESMF_Grid) :: grid real(ESMF_KIND_R4), allocatable :: valueListr4(:) real(ESMF_KIND_R8), allocatable :: valueListr8(:) @@ -3088,6 +3190,8 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) logical :: thereAreVerticals integer :: ch_dimid, timeiso_varid character(len=ESMF_MAXSTR) :: time_iso + integer :: wrt_mpi_comm + type(ESMF_VM) :: vm rc = ESMF_SUCCESS @@ -3132,6 +3236,26 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) status = ESMF_FILESTATUS_REPLACE endif + if ( tileFileName(1:7) == 'RESTART' ) then + + ! Write out restart files using write_restart_netcdf, then return from this subroutine + if (timeslice == 1) then + call ESMF_GridCompGet(comp, localPet=localPet, petCount=petCount, vm=vm, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_VMGet(vm=vm, mpiCommunicator=wrt_mpi_comm, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (petCount > 1) then + call write_restart_netcdf(wrtTileFB, trim(tileFileName), .true., wrt_mpi_comm, localPet, rc) + else + call write_restart_netcdf(wrtTileFB, trim(tileFileName), .false., wrt_mpi_comm, localPet, rc) + endif + + endif + return + endif + call ESMF_LogWrite("In ioCompRun() before writing to: "// & trim(tileFileName), ESMF_LOGMSG_INFO, rc=rc) @@ -3173,7 +3297,7 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) if (.not.isPresent) cycle ! field does not have the AttPack call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3", & name="ESMF:ungridded_dim_labels", isPresent=isPresent, & - itemCount=udimCount, rc=rc) + itemCount=udimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -3189,7 +3313,10 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) ! loop over all ungridded dimension labels do k=1, udimCount call write_out_ungridded_dim_atts(dimLabel=trim(udimList(k)), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! for restart files we store ungridded dimension labels in fields + call write_out_ungridded_dim_atts_from_field(field, dimLabel=trim(udimList(k)), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo deallocate(udimList) @@ -3197,7 +3324,7 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) deallocate(fieldList) if (thereAreVerticals) then ! see if the vertical_dim_labels attribute exists on the grid, and - ! if so access it and write out vecticals accordingly + ! if so access it and write out verticals accordingly call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name="vertical_dim_labels", isPresent=isPresent, & itemCount=udimCount, rc=rc) @@ -3259,15 +3386,12 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_enddef(ncid=ncid) - if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, varid, values=time) - if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ncerr = nf90_put_var(ncid, timeiso_varid, values=[trim(time_iso)]) - if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return ! loop over all the grid attributes that start with "time:", and @@ -3409,15 +3533,19 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc) character(len=*) :: dimLabel integer, intent(out) :: rc + logical :: isPresent + ! inquire if NetCDF file already contains this ungridded dimension ncerr = nf90_inq_varid(ncid, trim(dimLabel), varid=varid) if (ncerr == NF90_NOERR) return ! the variable does not exist in the NetCDF file yet -> add it ! access the undistributed dimension attribute on the grid call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & - name=trim(dimLabel), itemCount=valueCount, typekind=typekind, rc=rc) + name=trim(dimLabel), isPresent=isPresent, itemCount=valueCount, typekind=typekind, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (.not.isPresent) return ! nothing there to do + if( typekind == ESMF_TYPEKIND_R4 ) then allocate(valueListr4(valueCount)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & @@ -3430,6 +3558,8 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc) name=trim(dimLabel), valueList=valueListr8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + write(0,*) 'in write_out_ungridded_dim_atts: ERROR unknown typekind' endif ! now add it to the NetCDF file ncerr = nf90_redef(ncid=ncid) @@ -3540,6 +3670,147 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc) ncerr = nf90_enddef(ncid=ncid) if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return endif + end subroutine write_out_ungridded_dim_atts + + subroutine write_out_ungridded_dim_atts_from_field(field, dimLabel, rc) + + type(ESMF_Field),intent(in) :: field + character(len=*),intent(in) :: dimLabel + integer, intent(out) :: rc + + ! inquire if NetCDF file already contains this ungridded dimension + ncerr = nf90_inq_varid(ncid, trim(dimLabel), varid=varid) + if (ncerr == NF90_NOERR) return + ! the variable does not exist in the NetCDF file yet -> add it + ! access the undistributed dimension attribute on the grid + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", & + name=trim(dimLabel), itemCount=valueCount, typekind=typekind, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if( typekind == ESMF_TYPEKIND_R4 ) then + allocate(valueListr4(valueCount)) + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", & + name=trim(dimLabel), valueList=valueListr4, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + else if ( typekind == ESMF_TYPEKIND_R8) then + allocate(valueListr8(valueCount)) + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", & + name=trim(dimLabel), valueList=valueListr8, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + write(0,*) 'in write_out_ungridded_dim_atts: ERROR unknown typekind' + endif + ! now add it to the NetCDF file + ncerr = nf90_redef(ncid=ncid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + ncerr = nf90_inq_dimid(ncid, trim(dimLabel), dimid=dimid) + if (ncerr /= NF90_NOERR) then + ! dimension does not yet exist, and must be defined + ncerr = nf90_def_dim(ncid, trim(dimLabel), valueCount, dimid=dimid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + endif + if( typekind == ESMF_TYPEKIND_R4 ) then + ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_FLOAT, & + dimids=(/dimid/), varid=varid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + ncerr = nf90_enddef(ncid=ncid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + ncerr = nf90_put_var(ncid, varid, values=valueListr4) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + deallocate(valueListr4) + else if(typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_DOUBLE, & + dimids=(/dimid/), varid=varid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + ncerr = nf90_enddef(ncid=ncid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + ncerr = nf90_put_var(ncid, varid, values=valueListr8) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + deallocate(valueListr8) + endif + ! add attributes to this vertical variable + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", & + attnestflag=ESMF_ATTNEST_OFF, count=attCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (attCount>0) then + ncerr = nf90_redef(ncid=ncid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + endif + ! loop over all the attributes + do j=1, attCount + call ESMF_AttributeGet(field, convention="NetCDF", purpose="FV3-dim", & + attnestflag=ESMF_ATTNEST_OFF, attributeIndex=j, & + name=attName, typekind=typekind, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ! test for name starting with trim(dimLabel)":" + if (index(trim(attName), trim(dimLabel)//":") == 1) then + ind = len(trim(dimLabel)//":") + ! found a matching attributes + if (typekind == ESMF_TYPEKIND_CHARACTER) then + call ESMF_AttributeGet(field, & + convention="NetCDF", purpose="FV3-dim", & + name=trim(attName), value=valueS, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ncerr = nf90_put_att(ncid, varid, & + trim(attName(ind+1:len(attName))), values=valueS) + + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + else if (typekind == ESMF_TYPEKIND_I4) then + call ESMF_AttributeGet(field, & + convention="NetCDF", purpose="FV3-dim", & + name=trim(attName), value=valueI4, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ncerr = nf90_put_att(ncid, varid, & + trim(attName(ind+1:len(attName))), values=valueI4) + + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + else if (typekind == ESMF_TYPEKIND_R4) then + call ESMF_AttributeGet(field, & + convention="NetCDF", purpose="FV3-dim", & + name=trim(attName), value=valueR4, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ncerr = nf90_put_att(ncid, varid, & + trim(attName(ind+1:len(attName))), values=valueR4) + + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + else if (typekind == ESMF_TYPEKIND_R8) then + call ESMF_AttributeGet(field, & + convention="NetCDF", purpose="FV3-dim", & + name=trim(attName), value=valueR8, rc=rc) + + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + ncerr = nf90_put_att(ncid, varid, & + trim(attName(ind+1:len(attName))), values=valueR8) + + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + endif + endif + enddo + if (attCount>0) then + ncerr = nf90_enddef(ncid=ncid) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + endif end subroutine end subroutine ioCompRun @@ -3563,6 +3834,7 @@ subroutine ESMFproto_FieldMakeSingleTile(field, tile, tileField, petList, rc) ! The original field passed in remains valid. type(ESMF_TypeKind_Flag) :: typekind + type(ESMF_StaggerLoc) :: staggerloc type(ESMF_Index_Flag) :: indexflag type(ESMF_Grid) :: grid, tileGrid type(ESMF_Array) :: array @@ -3594,6 +3866,7 @@ subroutine ESMFproto_FieldMakeSingleTile(field, tile, tileField, petList, rc) ! access information from the incoming field call ESMF_FieldGet(field, array=array, typekind=typekind, & + staggerloc=staggerloc, & dimCount=fieldDimCount, name=fieldName, grid=grid, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -3723,6 +3996,7 @@ subroutine ESMFproto_FieldMakeSingleTile(field, tile, tileField, petList, rc) ! create the tile specific field from the array tileField = ESMF_FieldCreate(tileGrid, array=array, name=fieldName, & + ! staggerloc=staggerloc, & gridToFieldMap=gridToFieldMap, ungriddedLBound=ungriddedLBound, & ungriddedUBound=ungriddedUBound, rc=rc) @@ -4079,7 +4353,7 @@ end subroutine lambert ! !----------------------------------------------------------------------- ! - subroutine get_outfile(nfl, filename, outfile_name,noutfile) + subroutine get_outfile(nfl, filename, outfile_name, noutfile) integer, intent(in) :: nfl character(*), intent(in) :: filename(:,:) character(*), intent(inout) :: outfile_name(:) diff --git a/io/post_fv3.F90 b/io/post_fv3.F90 index c76431a9b..33f098ab9 100644 --- a/io/post_fv3.F90 +++ b/io/post_fv3.F90 @@ -247,12 +247,18 @@ subroutine post_getattr_fv3(wrt_int_state,grid_id) real(4), dimension(:), allocatable :: ak4,bk4 real(8), dimension(:), allocatable :: ak8,bk8 type(ESMF_FieldBundle) :: fldbundle + character(128) :: wrtFBName ! spval = 9.99e20 ! field bundle do nfb=1, wrt_int_state%FBcount fldbundle = wrt_int_state%wrtFB(nfb) + call ESMF_FieldBundleGet(fldbundle, name=wrtFBName, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__,file=__FILE__)) return + + if (wrtFBName(1:8) == 'restart_') cycle + ! set grid spec: ! if(mype==0) print*,'in post_getattr_lam,output_grid=',trim(output_grid(grid_id)),'nfb=',nfb ! if(mype==0) print*,'in post_getattr_lam, lon1=',lon1,lon2,lat1,lat2,dlon,dlat @@ -861,6 +867,12 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp) get_lsmsk: do ibdl=1, wrt_int_state%FBCount + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(ibdl), name=wrtFBName, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__,file=__FILE__)) return + + if (wrtFBName(1:8) == 'restart_') cycle + + call ESMF_AttributeGet(wrt_int_state%wrtFB(ibdl), convention="NetCDF", purpose="FV3", & name="grid_id", value=bundle_grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -944,6 +956,13 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp) ! if(mype==0) print *,'in setvar, read field, ibdl=',ibdl,'idim=', & ! ista,iend,'jdim=',jsta,jend + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(ibdl), grid=wrtGrid, & + fieldCount=ncount_field, name=wrtFBName,rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + + if (wrtFBName(1:8) == 'restart_') cycle + call ESMF_AttributeGet(wrt_int_state%wrtFB(ibdl), convention="NetCDF", purpose="FV3", & name="grid_id", value=bundle_grid_id, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -951,10 +970,6 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp) if (grid_id /= bundle_grid_id) cycle - call ESMF_FieldBundleGet(wrt_int_state%wrtFB(ibdl), grid=wrtGrid, & - fieldCount=ncount_field, name=wrtFBName,rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, file=__FILE__)) return ! bail out call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90 index c8df6acc0..3897ff43c 100644 --- a/module_fcst_grid_comp.F90 +++ b/module_fcst_grid_comp.F90 @@ -56,6 +56,10 @@ module module_fcst_grid_comp use data_override_mod, only: data_override_init use fv_nggps_diags_mod, only: fv_dyn_bundle_setup use fv3gfs_io_mod, only: fv_phys_bundle_setup + use fv3gfs_restart_io_mod, only: fv_phy_restart_bundle_setup, fv_sfc_restart_bundle_setup + use fv_ufs_restart_io_mod, only: fv_core_restart_bundle_setup, & + fv_srf_wnd_restart_bundle_setup, & + fv_tracer_restart_bundle_setup use fms2_io_mod, only: FmsNetcdfFile_t, open_file, close_file, variable_exists, read_data @@ -64,7 +68,8 @@ module module_fcst_grid_comp use module_fv3_io_def, only: num_pes_fcst, num_files, filename_base, & nbdlphys, iau_offset use module_fv3_config, only: dt_atmos, fcst_mpi_comm, fcst_ntasks, & - quilting, calendar, cpl_grid_id, & + quilting, quilting_restart, & + calendar, cpl_grid_id, & cplprint_flag use get_stochy_pattern_mod, only: write_stoch_restart_atm @@ -73,6 +78,7 @@ module module_fcst_grid_comp use module_cplfields, only: realizeConnectedCplFields use atmos_model_mod, only: setup_exportdata + use CCPP_data, only: GFS_control ! !----------------------------------------------------------------------- ! @@ -100,6 +106,8 @@ module module_fcst_grid_comp integer :: numTracers = 0 integer :: frestart(999) + + integer :: mype ! !----------------------------------------------------------------------- ! @@ -310,7 +318,7 @@ subroutine init_dyn_fb(nest, importState, exportState, clock, rc) type(ESMF_Grid) :: grid integer :: itemCount - character(len=ESMF_MAXSTR) :: itemNameList(1) + character(len=ESMF_MAXSTR) :: itemNameList(1), fb_name type(ESMF_FieldBundle) :: fb, fcstFB call ESMF_GridCompGet(nest, grid=grid, rc=rc) @@ -340,9 +348,23 @@ subroutine init_dyn_fb(nest, importState, exportState, clock, rc) call ESMF_StateAdd(exportState,(/fb/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call fv_dyn_bundle_setup(Atmos%axes, fb, grid, quilting=.true., rc=rc) + call ESMF_FieldBundleGet(fb, name=fb_name, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (fb_name(1:19) == 'restart_fv_core.res') then + call fv_core_restart_bundle_setup(fb, grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else if (fb_name(1:22) == 'restart_fv_srf_wnd.res') then + call fv_srf_wnd_restart_bundle_setup(fb, grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else if (fb_name(1:21) == 'restart_fv_tracer.res') then + call fv_tracer_restart_bundle_setup(fb, grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + call fv_dyn_bundle_setup(Atmos%axes, fb, grid, quilting=.true., rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + endif + end subroutine init_dyn_fb ! !----------------------------------------------------------------------- @@ -359,6 +381,7 @@ subroutine init_phys_fb(nest, importState, exportState, clock, rc) type(ESMF_Grid) :: grid integer :: itemCount, i character(len=ESMF_MAXSTR), allocatable :: itemNameList(:) + character(len=ESMF_MAXSTR) :: fb_name type(ESMF_FieldBundle), allocatable :: fbList(:) type(ESMF_FieldBundle) :: fcstFB @@ -387,9 +410,21 @@ subroutine init_phys_fb(nest, importState, exportState, clock, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo - call fv_phys_bundle_setup(Atmos%diag, Atmos%axes, fbList, grid, quilting=.true., nbdlphys=itemCount, rc=rc) + ! get the name of the first field bundle and based on that determine if it's a history or restart bundles + call ESMF_FieldBundleGet(fbList(1), name=fb_name, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (fb_name(1:16) == 'restart_phy_data') then + call fv_phy_restart_bundle_setup(fbList(1), grid, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + elseif (fb_name(1:16) == 'restart_sfc_data') then + call fv_sfc_restart_bundle_setup(fbList(1), grid, GFS_control, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else + call fv_phys_bundle_setup(Atmos%diag, Atmos%axes, fbList, grid, quilting=.true., nbdlphys=itemCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + endif + end subroutine init_phys_fb ! !----------------------------------------------------------------------- @@ -497,7 +532,6 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) integer,dimension(6) :: date, date_end ! integer :: initClock, unit, total_inttime - integer :: mype integer :: stat character(4) dateSY character(2) dateSM,dateSD,dateSH,dateSN,dateSS @@ -507,8 +541,10 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) character(256) :: gridfile character(8) :: bundle_grid - type(ESMF_FieldBundle),dimension(:), allocatable :: fieldbundle ! dynamics bundles - type(ESMF_FieldBundle),dimension(:,:), allocatable :: fieldbundlephys ! physics bundles + type(ESMF_FieldBundle),dimension(:), allocatable :: fieldbundle ! dynamics hystory bundles + type(ESMF_FieldBundle),dimension(:,:), allocatable :: fieldbundle_dyn_restart ! dynamics restart bundles + type(ESMF_FieldBundle),dimension(:,:), allocatable :: fieldbundlephys ! physics hystory bundles + type(ESMF_FieldBundle),dimension(:,:), allocatable :: fieldbundle_phy_restart ! physics restart bundles real(kind=8) :: mpi_wtime, timeis @@ -543,6 +579,7 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) integer,allocatable :: grid_number_on_all_pets(:) logical,allocatable :: is_moving_on_all_pets(:), is_moving(:) + character(len=7) :: nest_suffix type(FmsNetcdfFile_t) :: fileobj ! @@ -952,6 +989,9 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) nbdlphys = 2 allocate(fieldbundlephys(nbdlphys,ngrids)) + allocate(fieldbundle_dyn_restart(ngrids,3)) ! fv_core.res fv_srf_wnd.res fv_tracer.res + allocate(fieldbundle_phy_restart(ngrids,2)) ! phy_data sfc_data + do n=1,ngrids bundle_grid='' if (ngrids > 1 .and. n >= 2) then @@ -979,6 +1019,14 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) name="grid_id", value=n, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_AttributeAdd(fieldbundle(n), convention="NetCDF", purpose="FV3-nooutput", & + attrList=(/"frestart"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundle(n), convention="NetCDF", purpose="FV3-nooutput", & + name="frestart", valueList=frestart, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_StateAdd(tempState, (/fieldbundle(n)/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return @@ -1006,6 +1054,14 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) name="grid_id", value=n, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_AttributeAdd(fieldbundlephys(j,n), convention="NetCDF", purpose="FV3-nooutput", & + attrList=(/"frestart"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundlephys(j,n), convention="NetCDF", purpose="FV3-nooutput", & + name="frestart", valueList=frestart, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_StateAdd(tempState, (/fieldbundlephys(j,n)/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo @@ -1025,9 +1081,120 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) call ESMF_StateDestroy(tempState, noGarbage=.true., rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - enddo ! num_files + enddo ! num_files history + + if ( quilting_restart ) then + + do i=1,3 ! 3 dynamics restart bundles + + tempState = ESMF_StateCreate(rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (i == 1) then + name_FB = 'restart_fv_core.res' + elseif (i == 2) then + name_FB = 'restart_fv_srf_wnd.res' + elseif (i == 3) then + name_FB = 'restart_fv_tracer.res' + else + write(0,*)' unknown name_dynamics restart bundle ', i + ESMF_ERR_ABORT(101) + endif + + if (n > 1) then + write(nest_suffix,'(A5,I2.2)') '.nest', n + name_FB = trim(name_FB)//nest_suffix + endif + + fieldbundle_dyn_restart(n,i) = ESMF_FieldBundleCreate(name=trim(name_FB),rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeAdd(fieldbundle_dyn_restart(n,i), convention="NetCDF", purpose="FV3", & + attrList=(/"grid_id"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundle_dyn_restart(n,i), convention="NetCDF", purpose="FV3", & + name="grid_id", value=n, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeAdd(fieldbundle_dyn_restart(n,i), convention="NetCDF", purpose="FV3-nooutput", & + attrList=(/"frestart"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundle_dyn_restart(n,i), convention="NetCDF", purpose="FV3-nooutput", & + name="frestart", valueList=frestart, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_StateAdd(tempState, (/fieldbundle_dyn_restart(n,i)/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_GridCompInitialize(fcstGridComp(n), importState=tempState, & + exportState=exportState, phase=1, userrc=urc, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ESMF_LogFoundError(rcToCheck=urc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call ESMF_StateDestroy(tempState, noGarbage=.true., rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + enddo ! 3 dynamics restart bundles + + do i=1,2 ! 2 physics restart bundles + + tempState = ESMF_StateCreate(rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + if (i == 1) then + name_FB = 'restart_phy_data' + elseif (i == 2) then + name_FB = 'restart_sfc_data' + else + write(0,*)' unknown name_physics restart bundle ', i + ESMF_ERR_ABORT(101) + endif + + if (n > 1) then + write(nest_suffix,'(A5,I2.2)') '.nest', n + name_FB = trim(name_FB)//nest_suffix + endif + + fieldbundle_phy_restart(n,i) = ESMF_FieldBundleCreate(name=trim(name_FB),rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeAdd(fieldbundle_phy_restart(n,i), convention="NetCDF", purpose="FV3", & + attrList=(/"grid_id"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundle_phy_restart(n,i), convention="NetCDF", purpose="FV3", & + name="grid_id", value=n, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeAdd(fieldbundle_phy_restart(n,i), convention="NetCDF", purpose="FV3-nooutput", & + attrList=(/"frestart"/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fieldbundle_phy_restart(n,i), convention="NetCDF", purpose="FV3-nooutput", & + name="frestart", valueList=frestart, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_StateAdd(tempState, (/fieldbundle_phy_restart(n,i)/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_GridCompInitialize(fcstGridComp(n), importState=tempState, & + exportState=exportState, phase=2, userrc=urc, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ESMF_LogFoundError(rcToCheck=urc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + + call ESMF_StateDestroy(tempState, noGarbage=.true., rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + enddo ! 2 physics restart bundles + endif ! quilting_restart + enddo ! ngrids + ! total number of field bundles created is ngrids * (1(atm) + 2(phy) + 3(dyn_rest) +2(phy_rest) + if (mype == 0) write(*,*)'fcst_initialize: total number of field bundles: ', ngrids*(1+2+0+2) + !end qulting endif @@ -1058,7 +1225,6 @@ subroutine fcst_advertise(fcst_comp, importState, exportState, clock, rc) ! !*** local variables type(ESMF_VM) :: vm - integer :: mype integer :: n integer :: urc @@ -1070,7 +1236,7 @@ subroutine fcst_advertise(fcst_comp, importState, exportState, clock, rc) call ESMF_VMGetCurrent(vm=vm,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_VMGet(vm=vm, localPet=mype, rc=rc) + call ESMF_VMGet(vm=vm, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (mype == 0) write(*,*)'fcst_advertise, cpl_grid_id=',cpl_grid_id @@ -1100,7 +1266,6 @@ subroutine fcst_realize(fcst_comp, importState, exportState, clock, rc) ! !*** local variables type(ESMF_VM) :: vm - integer :: mype integer :: n integer :: urc @@ -1112,7 +1277,7 @@ subroutine fcst_realize(fcst_comp, importState, exportState, clock, rc) call ESMF_VMGetCurrent(vm=vm,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_VMGet(vm=vm, localPet=mype, rc=rc) + call ESMF_VMGet(vm=vm, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return if (mype == 0) write(*,*)'fcst_realize, cpl_grid_id=',cpl_grid_id @@ -1146,7 +1311,7 @@ subroutine fcst_run_phase_1(fcst_comp, importState, exportState,clock,rc) logical,save :: first=.true. integer,save :: dt_cap=0 type(ESMF_Time) :: currTime,stopTime - integer :: mype, seconds + integer :: seconds real(kind=8) :: mpi_wtime, tbeg1 ! !----------------------------------------------------------------------- @@ -1158,9 +1323,6 @@ subroutine fcst_run_phase_1(fcst_comp, importState, exportState,clock,rc) ! !----------------------------------------------------------------------- ! - call ESMF_GridCompGet(fcst_comp, localpet=mype, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call get_time(Atmos%Time - Atmos%Time_init, seconds) n_atmsteps = seconds/dt_atmos @@ -1214,7 +1376,7 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) ! !*** local variables ! - integer :: mype, date(6), seconds + integer :: date(6), seconds character(len=64) :: timestamp integer :: unit real(kind=8) :: mpi_wtime, tbeg1 @@ -1228,10 +1390,6 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) ! !----------------------------------------------------------------------- ! - call ESMF_GridCompGet(fcst_comp, localpet=mype, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return -! -!----------------------------------------------------------------------- ! *** call fcst integration subroutines call atmos_model_exchange_phase_2 (Atmos, rc=rc) @@ -1252,8 +1410,7 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc) !----- write restart file ------ if (mpp_pe() == mpp_root_pe())then - call get_date (Atmos%Time, date(1), date(2), date(3), & - date(4), date(5), date(6)) + call get_date (Atmos%Time, date(1), date(2), date(3), date(4), date(5), date(6)) open( newunit=unit, file='RESTART/'//trim(timestamp)//'.coupler.res' ) write( unit, '(i6,8x,a)' )calendar_type, & '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' @@ -1290,7 +1447,6 @@ subroutine fcst_finalize(fcst_comp, importState, exportState,clock,rc) ! !*** local variables ! - integer :: mype integer :: unit integer,dimension(6) :: date real(kind=8) :: mpi_wtime, tbeg1 @@ -1302,9 +1458,6 @@ subroutine fcst_finalize(fcst_comp, importState, exportState,clock,rc) tbeg1 = mpi_wtime() rc = ESMF_SUCCESS - call ESMF_GridCompGet(fcst_comp, localpet=mype, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call atmos_model_end (Atmos) call diag_manager_end (Atmos%Time) diff --git a/module_fv3_config.F90 b/module_fv3_config.F90 index a62800b34..e578e503e 100644 --- a/module_fv3_config.F90 +++ b/module_fv3_config.F90 @@ -19,7 +19,7 @@ module module_fv3_config ! integer :: cpl_grid_id logical :: cplprint_flag - logical :: quilting, output_1st_tstep_rst + logical :: quilting, quilting_restart, output_1st_tstep_rst ! real,dimension(:),allocatable :: output_fh character(esmf_maxstr),dimension(:),allocatable :: filename_base