From b8f7b0f33df6dc712f241dcd6818c9ec27134e26 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Wed, 7 Feb 2024 16:23:36 -0500 Subject: [PATCH 1/7] Update i/o routines for new surface coldstart file. (#698) * Update read of surface restart file to make the surface-specific records mandatory and the 'composite' records optional. Compute the 'composite' records from the surface-specific records. * Update logic to read either an 'old' style surface coldstart file or a new surface file from the new chgres_cube. The new files will have a global attribute saying this is a version 2 file. --- io/fv3atm_restart_io.F90 | 17 +- io/fv3atm_sfc_io.F90 | 398 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 393 insertions(+), 22 deletions(-) diff --git a/io/fv3atm_restart_io.F90 b/io/fv3atm_restart_io.F90 index d32be0586..a567f992f 100644 --- a/io/fv3atm_restart_io.F90 +++ b/io/fv3atm_restart_io.F90 @@ -14,7 +14,8 @@ module fv3atm_restart_io_mod register_axis, register_restart_field, & register_variable_attribute, register_field, & read_restart, write_restart, write_data, & - get_global_io_domain_indices, get_dimension_size + get_global_io_domain_indices, get_dimension_size, & + global_att_exists, get_global_attribute use mpp_domains_mod, only: domain2d use fv3atm_common_io, only: create_2d_field_and_add_to_bundle, & create_3d_field_and_add_to_bundle, copy_from_gfs_data, axis_type @@ -515,6 +516,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta !--- directory of the input files character(5) :: indir='INPUT' character(37) :: infile + character(2) :: file_ver !--- fms2_io file open logic logical :: amiopen logical :: override_frac_grid @@ -644,8 +646,19 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta amiopen=open_file(Sfc_restart, trim(infile), "read", domain=fv_domain, is_restart=.true., dont_add_res_to_filename=.true.) if( .not.amiopen ) call mpp_error(FATAL, 'Error opening file'//trim(infile)) + if (global_att_exists(Sfc_restart, "file_version")) then + call get_global_attribute(Sfc_restart, "file_version", file_ver) + if (file_ver == "V2") then + sfc%is_v2_file=.true. + endif + endif + if(sfc%allocate_arrays(Model, Atm_block, .true., warm_start)) then - call sfc%fill_2d_names(Model, warm_start) + if (sfc%is_v2_file) then + call sfc%fill_2d_names_v2(Model, warm_start) + else + call sfc%fill_2d_names(Model, warm_start) + endif call sfc%register_axes(Model, Sfc_restart, .true., warm_start) ! Tell CLM Lake to allocate data, and register its axes and fields diff --git a/io/fv3atm_sfc_io.F90 b/io/fv3atm_sfc_io.F90 index c0bfcf6d9..95957682a 100644 --- a/io/fv3atm_sfc_io.F90 +++ b/io/fv3atm_sfc_io.F90 @@ -21,7 +21,8 @@ module fv3atm_sfc_io private public :: Sfc_io_data_type - public :: Sfc_io_fill_2d_names, Sfc_io_fill_3d_names, Sfc_io_allocate_arrays, & + public :: Sfc_io_fill_2d_names, Sfc_io_fill_2d_names_v2, & + Sfc_io_fill_3d_names, Sfc_io_allocate_arrays, & Sfc_io_register_axes, Sfc_io_write_axes, Sfc_io_register_2d_fields, & Sfc_io_register_3d_fields, Sfc_io_copy_to_grid, Sfc_io_copy_from_grid, & Sfc_io_apply_safeguards, Sfc_io_transfer, Sfc_io_final @@ -49,6 +50,8 @@ module fv3atm_sfc_io ! The lsoil flag is only meaningful when reading:; logical, public :: is_lsoil = .false. + logical, public :: is_v2_file = .false. + ! SYNONYMS: Some nvar variables had two names in fv3atm_io.F90. They have ! only one name here. The "_s" is redundant because this file only has ! surface restart variables. @@ -79,6 +82,7 @@ module fv3atm_sfc_io procedure, public :: bundle_2d_fields => Sfc_io_bundle_2d_fields procedure, public :: bundle_3d_fields => Sfc_io_bundle_3d_fields procedure, public :: fill_2d_names => Sfc_io_fill_2d_names + procedure, public :: fill_2d_names_v2 => Sfc_io_fill_2d_names_v2 procedure, public :: fill_3d_names => Sfc_io_fill_3d_names procedure, public :: init_fields => Sfc_io_init_fields procedure, public :: transfer => Sfc_io_transfer @@ -558,6 +562,174 @@ subroutine Sfc_io_fill_2d_names(sfc,Model,warm_start) endif end subroutine Sfc_io_fill_2d_names + !>@ Fills the name2d array with all surface 2D field names. Updates nvar2m if needed. + !! This routine is for v2 coldstart files. + subroutine Sfc_io_fill_2d_names_v2(sfc,Model,warm_start) + implicit none + class(Sfc_io_data_type) :: sfc + type(GFS_control_type), intent(in) :: Model + 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 in file. Optional for cold starts. + nt=nt+1 ; sfc%name2(nt) = 'tg3' + nt=nt+1 ; sfc%name2(nt) = 'zorl' ! Optional for cold starts. + 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 in file. Optional for cold starts. + 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, unless indicated. + nt=nt+1 ; sfc%name2(nt) = 'scolor' + nt=nt+1 ; sfc%name2(nt) = 'sncovr' + nt=nt+1 ; sfc%name2(nt) = 'snodl' ! snowd on land portion of a cell. + ! Mandatory for cold starts. + nt=nt+1 ; sfc%name2(nt) = 'weasdl'! weasd on land portion of a cell. + ! Mandatory for cold starts. + 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 + ! Mandatory for cold starts. + 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 + ! Mandatory for cold starts. + 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. + ! Mandatory for cold starts. + nt=nt+1 ; sfc%name2(nt) = 'weasdi'! weasd on ice portion of a cell + ! Mandatory for cold starts. + + 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' + endif + + if(Model%cplwav) then + nt=nt+1 ; sfc%name2(nt) = 'zorlwav' !zorl from wave component + sfc%nvar2m = nt + endif + + 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' + endif + ! + ! 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 + + if (Model%lkm > 0 .and. Model%iopt_lake==Model%iopt_lake_flake) then + nt=nt+1 ; sfc%name2(nt) = 'T_snow' + nt=nt+1 ; sfc%name2(nt) = 'T_ice' + nt=nt+1 ; sfc%name2(nt) = 'h_ML' + nt=nt+1 ; sfc%name2(nt) = 't_ML' + nt=nt+1 ; sfc%name2(nt) = 't_mnw' + nt=nt+1 ; sfc%name2(nt) = 'h_talb' + nt=nt+1 ; sfc%name2(nt) = 't_talb' + nt=nt+1 ; sfc%name2(nt) = 't_bot1' + nt=nt+1 ; sfc%name2(nt) = 't_bot2' + nt=nt+1 ; sfc%name2(nt) = 'c_t' + endif + end subroutine Sfc_io_fill_2d_names_v2 + !>@ Registers 2D fields with FMS for reading or writing non-quilt restart files subroutine Sfc_io_register_2d_fields(sfc,Model,Sfc_restart,reading,warm_start) implicit none @@ -590,33 +762,62 @@ subroutine Sfc_io_register_2d_fields(sfc,Model,Sfc_restart,reading,warm_start) endif !--- register the 2D fields - do num = 1,sfc%nvar2m - var2_p => sfc%var2(:,:,num) - if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'tsfcl' .or. trim(sfc%name2(num)) == 'zorll' & - .or. trim(sfc%name2(num)) == 'zorli' .or. trim(sfc%name2(num)) == 'zorlwav' & - .or. trim(sfc%name2(num)) == 'snodl' .or. trim(sfc%name2(num)) == 'weasdl' & - .or. trim(sfc%name2(num)) == 'snodi' .or. trim(sfc%name2(num)) == 'weasdi' & - .or. trim(sfc%name2(num)) == 'tsfc' .or. trim(sfc%name2(num)) == 'zorlw' & + if (sfc%is_v2_file) then + do num = 1,sfc%nvar2m + var2_p => sfc%var2(:,:,num) + if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'zorll' & + .or. trim(sfc%name2(num)) == 'zorl' .or. trim(sfc%name2(num)) == 'zorlwav' & + .or. trim(sfc%name2(num)) == 'snwdph' .or. trim(sfc%name2(num)) == 'sheleg' & + .or. trim(sfc%name2(num)) == 'tsfc' & .or. trim(sfc%name2(num)) == 'albdirvis_lnd' .or. trim(sfc%name2(num)) == 'albdirnir_lnd' & .or. trim(sfc%name2(num)) == 'albdifvis_lnd' .or. trim(sfc%name2(num)) == 'albdifnir_lnd' & .or. trim(sfc%name2(num)) == 'albdirvis_ice' .or. trim(sfc%name2(num)) == 'albdirnir_ice' & .or. trim(sfc%name2(num)) == 'albdifvis_ice' .or. trim(sfc%name2(num)) == 'albdifnir_ice' & .or. trim(sfc%name2(num)) == 'emis_lnd' .or. trim(sfc%name2(num)) == 'emis_ice' & .or. trim(sfc%name2(num)) == 'sncovr_ice' .or. trim(sfc%name2(num)) == 'scolor') then - if(reading .and. sfc%is_lsoil) then - call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.) + if(reading .and. sfc%is_lsoil) then + call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.) + else + call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,& + & chunksizes=chunksizes2d, is_optional=.true.) + end if else - call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,& + if(reading .and. sfc%is_lsoil) then + call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/)) + else + call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d) + end if + endif + enddo + else + do num = 1,sfc%nvar2m + var2_p => sfc%var2(:,:,num) + if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'tsfcl' .or. trim(sfc%name2(num)) == 'zorll' & + .or. trim(sfc%name2(num)) == 'zorli' .or. trim(sfc%name2(num)) == 'zorlwav' & + .or. trim(sfc%name2(num)) == 'snodl' .or. trim(sfc%name2(num)) == 'weasdl' & + .or. trim(sfc%name2(num)) == 'snodi' .or. trim(sfc%name2(num)) == 'weasdi' & + .or. trim(sfc%name2(num)) == 'tsfc' .or. trim(sfc%name2(num)) == 'zorlw' & + .or. trim(sfc%name2(num)) == 'albdirvis_lnd' .or. trim(sfc%name2(num)) == 'albdirnir_lnd' & + .or. trim(sfc%name2(num)) == 'albdifvis_lnd' .or. trim(sfc%name2(num)) == 'albdifnir_lnd' & + .or. trim(sfc%name2(num)) == 'albdirvis_ice' .or. trim(sfc%name2(num)) == 'albdirnir_ice' & + .or. trim(sfc%name2(num)) == 'albdifvis_ice' .or. trim(sfc%name2(num)) == 'albdifnir_ice' & + .or. trim(sfc%name2(num)) == 'emis_lnd' .or. trim(sfc%name2(num)) == 'emis_ice' & + .or. trim(sfc%name2(num)) == 'sncovr_ice' .or. trim(sfc%name2(num)) == 'scolor') then + if(reading .and. sfc%is_lsoil) then + call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.) + else + call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,& & chunksizes=chunksizes2d, is_optional=.true.) - end if - else - if(reading .and. sfc%is_lsoil) then - call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/)) + end if else - call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d) - end if - endif - enddo + if(reading .and. sfc%is_lsoil) then + call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/)) + else + call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d) + end if + endif + enddo + endif if (Model%nstf_name(1) > 0) then mand = .false. @@ -702,7 +903,11 @@ subroutine Sfc_io_register_3d_fields(sfc,Model,Sfc_restart,reading,warm_start) !--- register the 3D fields var3_p => sfc%var3ice(:,:,:) - call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.true.) + if (sfc%is_v2_file) then + call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.false.) + else + call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.true.) + endif if(reading) then do num = 1,sfc%nvar3 @@ -1309,6 +1514,157 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop) enddo endif + if (sfc%is_v2_file) then + + if (sfc%var2(i,j,27) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snowd') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%fice(ix) > zero) then + Sfcprop(nb)%snowd(ix) = Sfcprop(nb)%snodi(ix) + elseif (Sfcprop(nb)%landfrac(ix) > zero) then + Sfcprop(nb)%snowd(ix) = Sfcprop(nb)%snodl(ix) + else + Sfcprop(nb)%snowd(ix) = zero + endif + enddo + enddo + endif + + if (sfc%var2(i,j,3) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing weasd') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%fice(ix) > zero) then + Sfcprop(nb)%weasd(ix) = Sfcprop(nb)%weasdi(ix) + elseif (Sfcprop(nb)%landfrac(ix) > zero) then + Sfcprop(nb)%weasd(ix) = Sfcprop(nb)%weasdl(ix) + else + Sfcprop(nb)%weasd(ix) = zero + endif + enddo + enddo + endif + +! Needed for first time step in radiation before Noah/NoahMP sets it from look up table. +! Just use a nominal value. + if (sfc%var2(i,j,39) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%landfrac(ix) > zero) then + Sfcprop(nb)%zorll(ix) = 25.0 + endif + enddo + enddo + endif + + if (sfc%var2(i,j,5) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorl') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%fice(ix) > zero) then + Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorli(ix) + elseif (Sfcprop(nb)%landfrac(ix) > zero) then + Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorll(ix) + else + Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlw(ix) + endif + enddo + enddo + endif + + if (sfc%var2(i,j,46) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing emis_ice') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%emis_ice(ix) = 0.96 + enddo + enddo + endif + + if (sfc%var2(i,j,47) < -9990.0_kind_phys .and. Model%lsm /= Model%lsm_ruc) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing sncovr_ice') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + ! Sfcprop(nb)%sncovr_ice(ix) = Sfcprop(nb)%sncovr(ix) + Sfcprop(nb)%sncovr_ice(ix) = zero + enddo + enddo + endif + + if (Model%use_cice_alb) then + if (sfc%var2(i,j,50) < -9990.0_kind_phys) then + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%oceanfrac(ix) > zero .and. & + Sfcprop(nb)%fice(ix) >= Model%min_seaice) then + Sfcprop(nb)%albdirvis_ice(ix) = 0.6_kind_phys + Sfcprop(nb)%albdifvis_ice(ix) = 0.6_kind_phys + Sfcprop(nb)%albdirnir_ice(ix) = 0.6_kind_phys + Sfcprop(nb)%albdifnir_ice(ix) = 0.6_kind_phys + endif + enddo + enddo + endif + + endif + + ! Fill in composite tsfc for coldstart runs - must happen after tsfcl is computed + compute_tsfc_for_coldstart: if (sfc%var2(i,j,36) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing composite tsfc') + if(Model%frac_grid) then ! 3-way composite + !$omp parallel do default(shared) private(nb, ix, tem, tem1) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) ! this may break restart reproducibility + tem1 = one - Sfcprop(nb)%landfrac(ix) + tem = tem1 * Sfcprop(nb)%fice(ix) ! tem = ice fraction wrt whole cell + Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix) * Sfcprop(nb)%landfrac(ix) & + + Sfcprop(nb)%tisfc(ix) * tem & + + Sfcprop(nb)%tsfco(ix) * (tem1-tem) + enddo + enddo + else + !$omp parallel do default(shared) private(nb, ix, tem) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%slmsk(ix) == 1) then + Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix) + if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc land ',nb,ix,Sfcprop(nb)%tsfcl(ix) + elseif(Sfcprop(nb)%fice(ix) > 0.0)then + Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tisfc(ix) + if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc ice ',nb,ix,Sfcprop(nb)%tisfc(ix) + else + Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix) + if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc water ',nb,ix,Sfcprop(nb)%tsfco(ix) + endif + enddo + enddo + endif + endif compute_tsfc_for_coldstart + + if (sfc%var2(i,j,sfc%nvar2m) < -9990.0_kind_phys) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorlwav') + !$omp parallel do default(shared) private(nb, ix) + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%zorlwav(ix) = Sfcprop(nb)%zorl(ix) !--- compute zorlwav from existing variables + enddo + enddo + endif + + + else ! old verion of coldstart file + + if (sfc%var2(i,j,34) < -9990.0_kind_phys) then if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snodl') !$omp parallel do default(shared) private(nb, ix, tem) @@ -1503,6 +1859,8 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop) enddo endif + endif ! check on which version of the surface file. + end subroutine Sfc_io_apply_safeguards !>@ destructor for Sfc_io_data_type From 9dec0faa6eafefc4b01f3cd117231619fb0670f4 Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Thu, 8 Feb 2024 16:23:33 -0500 Subject: [PATCH 2/7] update upp (#774) * Update UPP: fix the missing reflectivity bug --- upp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/upp b/upp index 78f369b01..945cb2cef 160000 --- a/upp +++ b/upp @@ -1 +1 @@ -Subproject commit 78f369b011ec41dca87a1260298d0228f2c4f38f +Subproject commit 945cb2cef5e8bd5949afd4f0fc35c4fb6e95a1bf From 28bfc365eb65364f1b4972a870837640cb7331a3 Mon Sep 17 00:00:00 2001 From: Dusan Jovic <48258889+DusanJovic-NOAA@users.noreply.github.com> Date: Fri, 9 Feb 2024 16:27:07 -0500 Subject: [PATCH 3/7] Fix out of bound errors in block_atmos_copy routines (#778) --- cpl/module_block_data.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cpl/module_block_data.F90 b/cpl/module_block_data.F90 index 7bd0a71b3..b5dc6cc80 100644 --- a/cpl/module_block_data.F90 +++ b/cpl/module_block_data.F90 @@ -365,7 +365,7 @@ subroutine block_array_copy_2d_r8_to_2d_r8(destin_ptr, source_arr, block, block_ jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j) = factor * source_arr(ib,jb) + destin_ptr(i,j) = factor * source_arr(i,j) enddo localrc = ESMF_SUCCESS end if @@ -440,7 +440,7 @@ subroutine block_array_copy_3d_r8_to_3d_r8(destin_ptr, source_arr, block, block_ jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j,k) = factor * source_arr(ib,jb,k) + destin_ptr(i,j,k) = factor * source_arr(i,j,k) enddo enddo localrc = ESMF_SUCCESS @@ -523,7 +523,7 @@ subroutine block_array_copy_3dslice_r8_to_3d_r8(destin_ptr, source_arr, slice, b jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j,k) = factor * source_arr(ib,jb,k,slice) + destin_ptr(i,j,k) = factor * source_arr(i,j,k,slice) enddo enddo localrc = ESMF_SUCCESS @@ -1004,7 +1004,7 @@ subroutine block_array_copy_2d_r4_to_2d_r8(destin_ptr, source_arr, block, block_ jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j) = factor * source_arr(ib,jb) + destin_ptr(i,j) = factor * source_arr(i,j) enddo localrc = ESMF_SUCCESS end if @@ -1079,7 +1079,7 @@ subroutine block_array_copy_3d_r4_to_3d_r8(destin_ptr, source_arr, block, block_ jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j,k) = factor * source_arr(ib,jb,k) + destin_ptr(i,j,k) = factor * source_arr(i,j,k) enddo enddo localrc = ESMF_SUCCESS @@ -1162,7 +1162,7 @@ subroutine block_array_copy_3dslice_r4_to_3d_r8(destin_ptr, source_arr, slice, b jb = block%index(block_index)%jj(ix) i = ib - block%isc + 1 j = jb - block%jsc + 1 - destin_ptr(i,j,k) = factor * source_arr(ib,jb,k,slice) + destin_ptr(i,j,k) = factor * source_arr(i,j,k,slice) enddo enddo localrc = ESMF_SUCCESS From 0fe9ba3020ad2f744ead7323bda414497aa33632 Mon Sep 17 00:00:00 2001 From: Dusan Jovic <48258889+DusanJovic-NOAA@users.noreply.github.com> Date: Wed, 14 Feb 2024 19:25:03 -0500 Subject: [PATCH 4/7] Add option to write 3d soilt, soilw and soill variables to a history file (#744) * Add 3d soil variables to history output * Flip only 3d variables with vertical dimension == levs (atm model levels) * Support 3d soil variables in inline post. * Update io/module_write_netcdf.F90 --- ccpp/driver/GFS_diagnostics.F90 | 140 +++++++++++++++++++++++--------- io/fv3atm_history_io.F90 | 116 ++++++++++++++++++-------- io/module_write_netcdf.F90 | 90 +++++++++++++++----- io/module_wrt_grid_comp.F90 | 18 ++++ io/post_fv3.F90 | 43 ++++++++++ 5 files changed, 314 insertions(+), 93 deletions(-) diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90 index 4c0e7cdf7..01d49593e 100644 --- a/ccpp/driver/GFS_diagnostics.F90 +++ b/ccpp/driver/GFS_diagnostics.F90 @@ -1,11 +1,11 @@ module GFS_diagnostics !----------------------------------------------------------------------- -! GFS_diagnostics_mod defines a data type and contains the routine +! GFS_diagnostics_mod defines a data type and contains the routine ! to populate said type with diagnostics from the GFS physics for ! use by the modeling system for output !----------------------------------------------------------------------- - + use machine, only: kind_phys !--- GFS_typedefs --- @@ -51,7 +51,7 @@ module GFS_diagnostics CONTAINS !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + ! Helper function for GFS_externaldiag_populate to handle the massive dtend(:,:,dtidx(:,:)) array subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit) implicit none @@ -62,7 +62,7 @@ subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit) integer, intent(inout) :: idx real(kind=kind_phys), pointer :: dtend(:,:,:) ! Assumption: dtend is null iff all(dtidx <= 1) character(len=*), intent(in), optional :: desc, unit - + integer :: idtend, nb idtend = Model%dtidx(itrac,iprocess) @@ -88,17 +88,17 @@ subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit) enddo endif end subroutine add_dtend - -!------------------------------------------------------------------------- + +!------------------------------------------------------------------------- !--- GFS_externaldiag_populate --- -!------------------------------------------------------------------------- -! creates and populates a data type with GFS physics diagnostic +!------------------------------------------------------------------------- +! creates and populates a data type with GFS physics diagnostic ! variables which is then handed off to the IPD for use by the model -! infrastructure layer to output as needed. The data type includes -! names, units, conversion factors, etc. There is no copying of data, -! but instead pointers are associated to the internal representation +! infrastructure layer to output as needed. The data type includes +! names, units, conversion factors, etc. There is no copying of data, +! but instead pointers are associated to the internal representation ! of each individual physics diagnostic. -!------------------------------------------------------------------------- +!------------------------------------------------------------------------- subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop, Coupling, & Grid, Tbd, Cldprop, Radtend, IntDiag, Init_parm) !---------------------------------------------------------------------------------------------! @@ -158,7 +158,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(:)%name = '' ExtDiag(:)%intpl_method = 'nearest_stod' - idx = 0 + idx = 0 idx = idx + 1 ExtDiag(idx)%axes = 2 @@ -949,7 +949,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ! !--- accumulated diagnostics --- do num = 1,NFXR - write (xtra,'(I2.2)') num + write (xtra,'(I2.2)') num idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'fluxr_'//trim(xtra) @@ -965,7 +965,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop !--- the next two appear to be appear to be coupling fields in gloopr !--- each has four elements !rab do num = 1,4 -!rab write (xtra,'(I1)') num +!rab write (xtra,'(I1)') num !rab idx = idx + 1 !rab ExtDiag(idx)%axes = 2 !rab ExtDiag(idx)%name = 'dswcmp_'//trim(xtra) @@ -978,7 +978,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop !rab enddo !rab !rab do num = 1,4 -!rab write (xtra,'(I1)') num +!rab write (xtra,'(I1)') num !rab idx = idx + 1 !rab ExtDiag(idx)%axes = 2 !rab ExtDiag(idx)%name = 'uswcmp_'//trim(xtra) @@ -1103,7 +1103,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snohfa(:) enddo - + if (Model%lsm == Model%lsm_noahmp) then idx = idx + 1 ExtDiag(idx)%axes = 2 @@ -1383,7 +1383,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%tedir(:) enddo - if (Model%lsm == Model%lsm_noahmp) then + if (Model%lsm == Model%lsm_noahmp) then idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'wa_acc' @@ -2197,7 +2197,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%smcref2(:) enddo - if (Model%lsm == Model%lsm_noahmp) then + if (Model%lsm == Model%lsm_noahmp) then idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'pahi' @@ -2468,7 +2468,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_pbl(:,:) enddo endif - + if (Model%do_spp) then idx = idx + 1 ExtDiag(idx)%axes = 3 @@ -2481,7 +2481,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_sfc(:,:) enddo endif - + if (Model%do_spp) then idx = idx + 1 ExtDiag(idx)%axes = 3 @@ -2494,7 +2494,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_mp(:,:) enddo endif - + if (Model%do_spp) then idx = idx + 1 ExtDiag(idx)%axes = 3 @@ -2507,7 +2507,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_gwd(:,:) enddo endif - + if (Model%do_spp) then idx = idx + 1 ExtDiag(idx)%axes = 3 @@ -2677,7 +2677,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%use_lake_model(:) enddo - + if(Model%iopt_lake==Model%iopt_lake_clm) then ! Populate the 3D arrays separately since the code is complicated: @@ -2704,7 +2704,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%lake_cannot_freeze(:) enddo - + idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'lake_t2m' @@ -2812,7 +2812,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%lake_ht(:) enddo - + endif endif @@ -2909,8 +2909,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%du3dt_pbl(:,:) enddo -! -! dv3dt_pbl +! +! dv3dt_pbl idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'dv3dt_pbl_ugwp' @@ -2921,8 +2921,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop do nb = 1,nblks ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%dv3dt_pbl(:,:) enddo -! -! dt3dt_pbl +! +! dt3dt_pbl idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'dt3dt_pbl_ugwp' @@ -2934,8 +2934,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%dt3dt_pbl(:,:) enddo ! -! uav_ugwp -! +! uav_ugwp +! idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'uav_ugwp' @@ -2947,8 +2947,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%uav_ugwp(:,:) enddo ! -! tav_ugwp -! +! tav_ugwp +! idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'tav_ugwp' @@ -2982,7 +2982,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%du3dt_ngw(:,:) enddo ! -! +! idx = idx + 1 ExtDiag(idx)%axes = 3 ExtDiag(idx)%name = 'du3dt_mtb' @@ -3454,7 +3454,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop endif enddo enddo - + if_qdiag3d: if(Model%qdiag3d) then idx = idx + 1 @@ -3499,7 +3499,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop !rab !rab do num = 1,5+Mdl_parms%pl_coeff -!rab write (xtra,'(I1)') num +!rab write (xtra,'(I1)') num !rab idx = idx + 1 !rab ExtDiag(idx)%axes = 3 !rab ExtDiag(idx)%name = 'dtend_'//trim(xtra) @@ -3877,7 +3877,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop idx = idx + 1 ExtDiag(idx)%axes = 2 ExtDiag(idx)%name = 'scolor' - ExtDiag(idx)%desc = 'soil color in integer 1-20' + ExtDiag(idx)%desc = 'soil color in integer 1-20' ExtDiag(idx)%unit = 'number' ExtDiag(idx)%mod_name = 'gfs_sfc' allocate (ExtDiag(idx)%data(nblks)) @@ -4203,6 +4203,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%sh2o(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soill' + ExtDiag(idx)%desc = 'liquid soil moisture' + ExtDiag(idx)%unit = 'm**3/m**3' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%sh2o(:,:) + enddo else do num = 1,Model%lsoil_lsm write (xtra,'(i1)') num @@ -4225,6 +4235,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%slc(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soill' + ExtDiag(idx)%desc = 'liquid soil moisture' + ExtDiag(idx)%unit = 'm**3/m**3' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%slc(:,:) + enddo endif if (Model%lsm == Model%lsm_ruc) then @@ -4241,6 +4261,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%smois(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soilw' + ExtDiag(idx)%desc = 'volumetric soil moisture' + ExtDiag(idx)%unit = 'fraction' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%smois(:,:) + enddo else do num = 1,Model%lsoil_lsm write (xtra,'(i1)') num @@ -4255,6 +4285,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%smc(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soilw' + ExtDiag(idx)%desc = 'volumetric soil moisture' + ExtDiag(idx)%unit = 'fraction' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%smc(:,:) + enddo endif if (Model%lsm == Model%lsm_ruc) then @@ -4271,6 +4311,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%tslb(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soilt' + ExtDiag(idx)%desc = 'soil temperature' + ExtDiag(idx)%unit = 'K' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%tslb(:,:) + enddo else do num = 1,Model%lsoil_lsm write (xtra,'(i1)') num @@ -4285,6 +4335,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%stc(:,num) enddo enddo + idx = idx + 1 + ExtDiag(idx)%axes = 3 + ExtDiag(idx)%name = 'soilt' + ExtDiag(idx)%desc = 'soil temperature' + ExtDiag(idx)%unit = 'K' + ExtDiag(idx)%mod_name = 'gfs_sfc' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%stc(:,:) + enddo endif !--------------------------nsst variables @@ -5396,7 +5456,7 @@ subroutine clm_lake_externaldiag_populate(ExtDiag, Model, Sfcprop, idx, cn_one, character(:), allocatable :: fullname integer :: nk, idx0, iblk - + do iblk=1,nblks call link_all_levels(Sfcprop(iblk)%lake_snow_z3d, 'lake_snow_z3d', 'lake snow level depth', 'm') enddo @@ -5524,6 +5584,6 @@ function soil_layer_depth(lsm, lsm_ruc, lsm_noah, layer) result(layer_depth) ! end function soil_layer_depth -!------------------------------------------------------------------------- +!------------------------------------------------------------------------- end module GFS_diagnostics diff --git a/io/fv3atm_history_io.F90 b/io/fv3atm_history_io.F90 index 7c73fe296..d7ee4e808 100644 --- a/io/fv3atm_history_io.F90 +++ b/io/fv3atm_history_io.F90 @@ -49,9 +49,10 @@ module fv3atm_history_io_mod type history_type integer :: tot_diag_idx = 0 - integer :: isco=0,ieco=0,jsco=0,jeco=0,levo=0,num_axes_phys=0 - integer :: fhzero=0, ncld=0, nsoil=0, imp_physics=0, landsfcmdl=0 + integer :: isco=0,ieco=0,jsco=0,jeco=0,num_axes_phys=0 + integer :: fhzero=0, ncld=0, nsoil=0, nsoil_lsm=0, imp_physics=0, landsfcmdl=0 real(4) :: dtp=0 + integer,dimension(:), pointer :: levo => null() integer,dimension(:), pointer :: nstt => null() integer,dimension(:), pointer :: nstt_vctbl => null() integer,dimension(:), pointer :: all_axes => null() @@ -182,11 +183,11 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat, hist%ieco = Atm_block%iec hist%jsco = Atm_block%jsc hist%jeco = Atm_block%jec - hist%levo = model%levs hist%fhzero = nint(Model%fhzero) ! hist%ncld = Model%ncld hist%ncld = Model%imp_physics hist%nsoil = Model%lsoil + hist%nsoil_lsm = Model%lsoil_lsm hist%dtp = Model%dtp hist%imp_physics = Model%imp_physics hist%landsfcmdl = Model%lsm @@ -208,7 +209,9 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat, call mpp_error(fatal, 'fv3atm_io::fv3atm_diag_register - need to increase parameter DIAG_SIZE') endif + allocate(hist%levo(hist%tot_diag_idx)) allocate(hist%nstt(hist%tot_diag_idx), hist%nstt_vctbl(hist%tot_diag_idx)) + hist%levo = 0 hist%nstt = 0 hist%nstt_vctbl = 0 nrgst_bl = 0 @@ -224,6 +227,7 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat, trim(Diag(idx)%unit), missing_value=real(missing_value)) if(Diag(idx)%id > 0) then if (Diag(idx)%axes == 2) then + hist%levo(idx) = 1 if( index(trim(Diag(idx)%intpl_method),'bilinear') > 0 ) then nrgst_bl = nrgst_bl + 1 hist%nstt(idx) = nrgst_bl @@ -239,17 +243,18 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat, endif endif else if (diag(idx)%axes == 3) then + hist%levo(idx) = size(Diag(idx)%data(1)%var3, dim=2) if( index(trim(diag(idx)%intpl_method),'bilinear') > 0 ) then hist%nstt(idx) = nrgst_bl + 1 - nrgst_bl = nrgst_bl + hist%levo + nrgst_bl = nrgst_bl + hist%levo(idx) else if (trim(diag(idx)%intpl_method) == 'nearest_stod' ) then hist%nstt(idx) = nrgst_nb + 1 - nrgst_nb = nrgst_nb + hist%levo + nrgst_nb = nrgst_nb + hist%levo(idx) endif if(trim(diag(idx)%intpl_method) == 'vector_bilinear') then if(diag(idx)%name(1:1) == 'v' .or. diag(idx)%name(1:1) == 'V') then hist%nstt_vctbl(idx) = nrgst_vctbl + 1 - nrgst_vctbl = nrgst_vctbl + hist%levo + nrgst_vctbl = nrgst_vctbl + hist%levo(idx) ! print *,'in phy_setup, vector_bilinear, name=', trim(diag(idx)%name),' nstt_vctbl=', hist%nstt_vctbl(idx), 'idx=',idx endif endif @@ -289,14 +294,14 @@ subroutine history_type_output(hist, time, diag, atm_block, nx, ny, levs, ntcw, real(kind=kind_phys), intent(in) :: time_radsw real(kind=kind_phys), intent(in) :: time_radlw !--- local variables - integer :: i, j, k, idx, nb, ix, ii, jj + integer :: i, j, k, idx, nb, ix, ii, jj, levo_3d character(len=2) :: xtra #ifdef CCPP_32BIT real, dimension(nx,ny) :: var2 - real, dimension(nx,ny,levs) :: var3 + real, dimension(:,:,:), allocatable :: var3 #else real(kind=kind_phys), dimension(nx,ny) :: var2 - real(kind=kind_phys), dimension(nx,ny,levs) :: var3 + real(kind=kind_phys), dimension(:,:,:), allocatable :: var3 #endif real(kind=kind_phys) :: rtime_int, rtime_intfull, lcnvfac real(kind=kind_phys) :: rtime_radsw, rtime_radlw @@ -444,31 +449,37 @@ subroutine history_type_output(hist, time, diag, atm_block, nx, ny, levs, ntcw, enddo endif if_mask endif int_or_real - ! used=send_data(Diag(idx)%id, var2, Time) - ! print *,'in phys, after store_data, idx=',idx,' var=', trim(Diag(idx)%name) + call hist%store_data(Diag(idx)%id, var2, Time, idx, Diag(idx)%intpl_method, Diag(idx)%name) - ! if(trim(Diag(idx)%name) == 'totprcp_ave' ) print *,'in gfs_io, totprcp=',Diag(idx)%data(1)%var2(1:3), & - ! ' lcnvfac=', lcnvfac + elseif (Diag(idx)%axes == 3) then !--- !--- skipping other 3D variables with the following else statement !--- - ! if(mpp_pe()==mpp_root_pe())print *,'in,fv3atm_io. 3D fields, idx=',idx,'varname=',trim(diag(idx)%name), & - ! 'lcnvfac=',lcnvfac, 'hist%levo=',hist%levo,'nx=',nx,'ny=',ny - do k=1, hist%levo + + levo_3d = hist%levo(idx) + allocate(var3(nx,ny,levo_3d)) + + do k=1, levo_3d do j = 1, ny jj = j + Atm_block%jsc -1 do i = 1, nx ii = i + Atm_block%isc -1 nb = Atm_block%blkno(ii,jj) ix = Atm_block%ixp(ii,jj) - ! if(mpp_pe()==mpp_root_pe())print *,'in,fv3atm_io,sze(Diag(idx)%data(nb)%var3)=', & - ! size(Diag(idx)%data(nb)%var3,1),size(Diag(idx)%data(nb)%var3,2) - var3(i,j,k) = Diag(idx)%data(nb)%var3(ix,hist%levo-k+1)*lcnvfac + ! flip only 3d variables with vertical dimension == levs (atm model levels) + if (levo_3d == levs) then + var3(i,j,k) = Diag(idx)%data(nb)%var3(ix,levo_3d-k+1)*lcnvfac + else + var3(i,j,k) = Diag(idx)%data(nb)%var3(ix, k )*lcnvfac + endif enddo enddo enddo + call hist%store_data3D(Diag(idx)%id, var3, Time, idx, Diag(idx)%intpl_method, Diag(idx)%name) + deallocate(var3) + endif if_2d endif has_id end do history_loop @@ -574,7 +585,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl #ifdef CCPP_32BIT real, intent(in) :: work(:,:,:) #else - real(kind=kind_phys), intent(in) :: work(hist%ieco-hist%isco+1,hist%jeco-hist%jsco+1,hist%levo) + real(kind=kind_phys), intent(in) :: work(hist%ieco-hist%isco+1,hist%jeco-hist%jsco+1,hist%levo(idx)) #endif type(time_type), intent(in) :: Time character(*), intent(in) :: intpl_method @@ -584,11 +595,12 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl integer k,j,i,nv,i1,j1 logical used ! + write(0,*) ' history_type_store_data3D kinds ', kind_phys, kind(work), lbound(work), ubound(work), size(work) if( id > 0 ) then if( hist%use_wrtgridcomp_output ) then if( trim(intpl_method) == 'bilinear') then !$omp parallel do default(shared) private(i,j,k) - do k= 1,hist%levo + do k= 1,hist%levo(idx) do j= hist%jsco,hist%jeco do i= hist%isco,hist%ieco hist%buffer_phys_bl(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k) @@ -597,7 +609,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl enddo else if(trim(intpl_method) == 'nearest_stod') then !$omp parallel do default(shared) private(i,j,k) - do k= 1,hist%levo + do k= 1,hist%levo(idx) do j= hist%jsco,hist%jeco do i= hist%isco,hist%ieco hist%buffer_phys_nb(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k) @@ -607,7 +619,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl else if(trim(intpl_method) == 'vector_bilinear') then !first save the data !$omp parallel do default(shared) private(i,j,k) - do k= 1,hist%levo + do k= 1,hist%levo(idx) do j= hist%jsco,hist%jeco do i= hist%isco,hist%ieco hist%buffer_phys_bl(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k) @@ -615,9 +627,9 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl enddo enddo if( fldname(1:1) == 'u' .or. fldname(1:1) == 'U') then - if(.not.associated(hist%uwork3d)) allocate(hist%uwork3d(hist%isco:hist%ieco,hist%jsco:hist%jeco,hist%levo)) + if(.not.associated(hist%uwork3d)) allocate(hist%uwork3d(hist%isco:hist%ieco,hist%jsco:hist%jeco,hist%levo(idx))) !$omp parallel do default(shared) private(i,j,k) - do k= 1, hist%levo + do k= 1, hist%levo(idx) do j= hist%jsco,hist%jeco do i= hist%isco,hist%ieco hist%uwork3d(i,j,k) = work(i-hist%isco+1,j-hist%jsco+1,k) @@ -642,7 +654,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl enddo enddo !$omp parallel do default(shared) private(i,j,k,nv,i1,j1) - do k= 1, hist%levo + do k= 1, hist%levo(idx) nv = hist%nstt_vctbl(idx)+k-1 do j= hist%jsco,hist%jeco j1 = j-hist%jsco+1 @@ -895,6 +907,36 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q deallocate(axis_data) enddo endif + + ! add zsoil axis + call ESMF_AttributeAdd(fcst_grid, convention="NetCDF", purpose="FV3", & + attrList=(/"zsoil ", & + "zsoil:long_name ", & + "zsoil:units ", & + "zsoil:cartesian_axis", & + "zsoil:positive "/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", & + name="zsoil", valueList=(/ (i, i=1,hist%nsoil_lsm) /), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", & + name="zsoil:long_name", value="soil level", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", & + name="zsoil:units", value="level", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", & + name="zsoil:cartesian_axis", value="Z", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", & + name="zsoil:positive", value="down", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + ! print *,'in setup fieldbundle_phys, hist%num_axes_phys=',hist%num_axes_phys,'hist%tot_diag_idx=',hist%tot_diag_idx, & ! 'nbdlphys=',nbdlphys ! @@ -913,7 +955,7 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q !add origin field call hist%add_field_to_phybundle(trim(output_name),trim(Diag(idx)%desc),trim(Diag(idx)%unit), "time: point", & - axes(1:Diag(idx)%axes), fcst_grid, hist%nstt(idx), phys_bundle(ibdl), outputfile(ibdl), & + axes(1:Diag(idx)%axes), fcst_grid, hist%nstt(idx), hist%levo(idx), phys_bundle(ibdl), outputfile(ibdl), & bdl_intplmethod(ibdl), rcd=rc) ! if( mpp_pe() == mpp_root_pe()) print *,'phys, add field,',trim(Diag(idx)%name),'idx=',idx,'ibdl=',ibdl ! @@ -923,7 +965,7 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q output_name = 'wind'//trim(output_name)//'vector' outputfile1 = 'none' call hist%add_field_to_phybundle(trim(output_name),trim(Diag(idx)%desc),trim(Diag(idx)%unit), "time: point", & - axes(1:Diag(idx)%axes), fcst_grid, hist%nstt_vctbl(idx),phys_bundle(ibdl), outputfile1, & + axes(1:Diag(idx)%axes), fcst_grid, hist%nstt_vctbl(idx), hist%levo(idx), phys_bundle(ibdl), outputfile1, & bdl_intplmethod(ibdl),l2dvector=l2dvector, rcd=rc) ! if( mpp_pe() == mpp_root_pe()) print *,'in phys, add vector field,',trim(Diag(idx)%name),' idx=',idx,' ibdl=',ibdl endif @@ -951,7 +993,7 @@ end subroutine history_type_bundle_setup !! It sets attributes for and logs information about a single ESMF field. Do not call this subroutine directly. !! Call fv_phys_bundle_setup instead. subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cell_methods, axes,phys_grid, & - kstt,phys_bundle,output_file,intpl_method,range,l2dvector,rcd) + kstt,levo,phys_bundle,output_file,intpl_method,range,l2dvector,rcd) ! use esmf ! @@ -961,7 +1003,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel character(*), intent(in) :: output_file, intpl_method integer, intent(in) :: axes(:) type(esmf_grid), intent(in) :: phys_grid - integer, intent(in) :: kstt + integer, intent(in) :: kstt, levo type(esmf_fieldbundle),intent(inout) :: phys_bundle real, intent(in), optional :: range(2) logical, intent(in), optional :: l2dvector @@ -1027,7 +1069,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) else if(size(axes) == 3) then - temp_r3d => hist%buffer_phys_nb(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+hist%levo-1) + temp_r3d => hist%buffer_phys_nb(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+levo-1) field = ESMF_FieldCreate(phys_grid, temp_r3d, datacopyflag=copyflag, & name=var_name, indexFlag=ESMF_INDEX_DELOCAL, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -1043,7 +1085,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) else if(size(axes) == 3) then - temp_r3d => hist%buffer_phys_bl(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+hist%levo-1) + temp_r3d => hist%buffer_phys_bl(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+levo-1) field = ESMF_FieldCreate(phys_grid, temp_r3d, datacopyflag=copyflag, & name=var_name, indexFlag=ESMF_INDEX_DELOCAL, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -1128,8 +1170,14 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", & attrList=(/"ESMF:ungridded_dim_labels"/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & - name="ESMF:ungridded_dim_labels", valueList=(/trim(hist%axis_name(idx))/), rc=rc) + + if (levo == hist%nsoil_lsm) then + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & + name="ESMF:ungridded_dim_labels", valueList=(/"zsoil"/), rc=rc) + else + call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & + name="ESMF:ungridded_dim_labels", valueList=(/trim(hist%axis_name(idx))/), rc=rc) + endif if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return endif enddo diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 082b38839..d9d8ff970 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -57,7 +57,7 @@ subroutine write_netcdf(wrtfb, filename, & ! !** local vars integer :: i,j,t, istart,iend,jstart,jend - integer :: im, jm, lm + integer :: im, jm, lm, lsoil integer, dimension(:), allocatable :: fldlev @@ -95,9 +95,10 @@ subroutine write_netcdf(wrtfb, filename, & integer :: ncerr,ierr integer :: ncid integer :: oldMode - integer :: im_dimid, jm_dimid, tile_dimid, pfull_dimid, phalf_dimid, time_dimid, ch_dimid + integer :: dim_len + integer :: im_dimid, jm_dimid, tile_dimid, pfull_dimid, phalf_dimid, time_dimid, ch_dimid, lsoil_dimid integer :: im_varid, jm_varid, tile_varid, lon_varid, lat_varid, timeiso_varid - integer, dimension(:), allocatable :: dimids_2d, dimids_3d, dimids, chunksizes + integer, dimension(:), allocatable :: dimids_2d, dimids_3d, dimids_soil, dimids, chunksizes integer, dimension(:), allocatable :: varids integer :: xtype integer :: quant_mode @@ -204,17 +205,24 @@ subroutine write_netcdf(wrtfb, filename, & end do lm = maxval(fldlev(:)) + call get_dimlen_if_exists(ncid, "zsoil", wrtgrid, lsoil, rc) + if (lsoil > 0 .and. (.not. any(fldlev(:) == lsoil))) then + lsoil = 0 + end if + if (lsoil > 0 .and. (.not. any(fldlev(:) > 1 .and. fldlev(:) /= lsoil))) then + lm = 1 + end if ! for serial output allocate 'global' arrays if (.not. par) then - allocate(array_r4(im,jm)) + ! allocate(array_r4(im,jm)) allocate(array_r8(im,jm)) - allocate(array_r4_3d(im,jm,lm)) + ! allocate(array_r4_3d(im,jm,lm)) allocate(array_r8_3d(im,jm,lm)) if (is_cubed_sphere) then - allocate(array_r4_cube(im,jm,tileCount)) + ! allocate(array_r4_cube(im,jm,tileCount)) allocate(array_r8_cube(im,jm,tileCount)) - allocate(array_r4_3d_cube(im,jm,lm,tileCount)) + ! allocate(array_r4_3d_cube(im,jm,lm,tileCount)) allocate(array_r8_3d_cube(im,jm,lm,tileCount)) end if end if @@ -243,6 +251,9 @@ subroutine write_netcdf(wrtfb, filename, & call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, mype, rc) call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, mype, rc) end if + if (lsoil > 1) then + call add_dim(ncid, "zsoil", lsoil_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 @@ -318,14 +329,16 @@ subroutine write_netcdf(wrtfb, filename, & ! define variables (fields) if (is_cubed_sphere) then allocate(dimids_2d(4)) - allocate(dimids_3d(5)) - dimids_2d = [im_dimid,jm_dimid, tile_dimid,time_dimid] - if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid,tile_dimid,time_dimid] + allocate(dimids_3d(5), dimids_soil(5)) + dimids_2d = [im_dimid,jm_dimid, tile_dimid,time_dimid] + if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid,tile_dimid,time_dimid] + if (lsoil > 1) dimids_soil = [im_dimid,jm_dimid,lsoil_dimid,tile_dimid,time_dimid] else allocate(dimids_2d(3)) - allocate(dimids_3d(4)) - dimids_2d = [im_dimid,jm_dimid, time_dimid] - if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid, time_dimid] + allocate(dimids_3d(4), dimids_soil(4)) + dimids_2d = [im_dimid,jm_dimid, time_dimid] + if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid, time_dimid] + if (lsoil > 1) dimids_soil = [im_dimid,jm_dimid,lsoil_dimid, time_dimid] end if do i=1, fieldCount @@ -336,7 +349,11 @@ subroutine write_netcdf(wrtfb, filename, & if (rank == 2) then dimids = dimids_2d else if (rank == 3) then - dimids = dimids_3d + if (fldlev(i) == lsoil) then + dimids = dimids_soil + else if (fldlev(i) == lm) then + dimids = dimids_3d + endif else if (mype==0) write(0,*)'Unsupported rank ', rank call ESMF_Finalize(endflag=ESMF_END_ABORT) @@ -368,7 +385,7 @@ subroutine write_netcdf(wrtfb, filename, & if (is_cubed_sphere) then chunksizes = [im, jm, lm, tileCount, 1] else - chunksizes = [ichunk3d(grid_id), jchunk3d(grid_id), kchunk3d(grid_id), 1] + chunksizes = [ichunk3d(grid_id), jchunk3d(grid_id), fldlev(i), 1] end if ncerr = nf90_def_var_chunking(ncid, varids(i), NF90_CHUNKED, chunksizes) ; NC_ERR_STOP(ncerr) end if @@ -610,17 +627,21 @@ subroutine write_netcdf(wrtfb, filename, & else if (is_cubed_sphere) then call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + allocate(array_r4_cube(im,jm,tileCount)) do t=1,tileCount call ESMF_ArrayGather(array, array_r4_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) end do if (do_io) then ncerr = nf90_put_var(ncid, varids(i), values=array_r4_cube, start=start_idx); NC_ERR_STOP(ncerr) end if + deallocate(array_r4_cube) 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 end if else if (typekind == ESMF_TYPEKIND_R8) then @@ -663,17 +684,21 @@ subroutine write_netcdf(wrtfb, filename, & else if (is_cubed_sphere) then call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + allocate(array_r4_3d_cube(im,jm,fldlev(i),tileCount)) do t=1,tileCount call ESMF_ArrayGather(array, array_r4_3d_cube(:,:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) end do if (mype==0) then ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d_cube, start=start_idx); NC_ERR_STOP(ncerr) end if + deallocate(array_r4_3d_cube) else + allocate(array_r4_3d(im,jm,fldlev(i))) 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 end if else if (typekind == ESMF_TYPEKIND_R8) then @@ -708,14 +733,14 @@ subroutine write_netcdf(wrtfb, filename, & end do ! end fieldCount if (.not. par) then - deallocate(array_r4) + ! deallocate(array_r4) deallocate(array_r8) - deallocate(array_r4_3d) + ! deallocate(array_r4_3d) deallocate(array_r8_3d) if (is_cubed_sphere) then - deallocate(array_r4_cube) + ! deallocate(array_r4_cube) deallocate(array_r8_cube) - deallocate(array_r4_3d_cube) + ! deallocate(array_r4_3d_cube) deallocate(array_r8_3d_cube) end if end if @@ -723,6 +748,7 @@ subroutine write_netcdf(wrtfb, filename, & if (do_io) then deallocate(dimids_2d) deallocate(dimids_3d) + deallocate(dimids_soil) end if deallocate(fcstField) @@ -883,6 +909,22 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) end subroutine get_grid_attr +!---------------------------------------------------------------------------------------- + subroutine get_dimlen_if_exists(ncid, dim_name, grid, dim_len, rc) + + integer, intent(in) :: ncid + character(len=*), intent(in) :: dim_name + type(ESMF_Grid), intent(in) :: grid + integer, intent(out) :: dim_len + integer, intent(out) :: rc + + dim_len = 0 + + call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & + name=dim_name, itemCount=dim_len, rc=rc); ESMF_ERR_RETURN(rc) + + end subroutine get_dimlen_if_exists + !> Add a dimension. !> !> @param[in] ncid NetCDF file ID. @@ -907,6 +949,7 @@ subroutine add_dim(ncid, dim_name, dimid, grid, mype, rc) integer :: ncerr type(ESMF_TypeKind_Flag) :: typekind + real(ESMF_KIND_I4), allocatable :: valueListI4(:) real(ESMF_KIND_R4), allocatable :: valueListR4(:) real(ESMF_KIND_R8), allocatable :: valueListR8(:) ! @@ -944,6 +987,15 @@ subroutine add_dim(ncid, dim_name, dimid, grid, mype, rc) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR4) + else if (typekind==ESMF_TYPEKIND_I4) then + ncerr = nf90_def_var(ncid, dim_name, NF90_INT4, dimids=[dimid], varid=dim_varid); NC_ERR_STOP(ncerr) + allocate(valueListI4(n)) + call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & + name=trim(dim_name), valueList=valueListI4, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, dim_varid, values=valueListI4); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + deallocate(valueListI4) else 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) diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index b7e93e28f..9a73f3d43 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -3363,6 +3363,7 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc) integer :: localPet, petCount, i, j, k, ind type(ESMF_Grid) :: grid + real(ESMF_KIND_I4), allocatable :: valueListi4(:) real(ESMF_KIND_R4), allocatable :: valueListr4(:) real(ESMF_KIND_R8), allocatable :: valueListr8(:) integer :: valueCount, fieldCount, udimCount @@ -3747,6 +3748,12 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dimLabel), valueList=valueListr8, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + else if ( typekind == ESMF_TYPEKIND_I4) then + allocate(valueListi4(valueCount)) + call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & + name=trim(dimLabel), valueList=valueListi4, 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' @@ -3784,6 +3791,17 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc) 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) + else if(typekind == ESMF_TYPEKIND_I4) then + ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_INT4, & + 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=valueListi4) + if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return + deallocate(valueListi4) endif ! add attributes to this vertical variable call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & diff --git a/io/post_fv3.F90 b/io/post_fv3.F90 index 97962bdd9..30eae6ba4 100644 --- a/io/post_fv3.F90 +++ b/io/post_fv3.F90 @@ -3562,6 +3562,49 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp) enddo endif + ! soilt + if(trim(fieldname)=='soilt') then + !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,stc,arrayr43d,sm,sice,fillvalue,spval) + do l=1,nsoil + do j=jsta,jend + do i=ista, iend + stc(i,j,l) = arrayr43d(i,j,l) + if( abs(arrayr43d(i,j,l)-fillValue) < small) stc(i,j,l) = spval + !mask open water areas, combine with sea ice tmp + if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,l) = spval + enddo + enddo + enddo + endif + + ! soilw + if(trim(fieldname)=='soilw') then + !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,smc,arrayr43d,sm,fillvalue,spval) + do l=1,nsoil + do j=jsta,jend + do i=ista, iend + smc(i,j,l) = arrayr43d(i,j,l) + if( abs(arrayr43d(i,j,l)-fillValue) < small) smc(i,j,l) = spval + if (sm(i,j) /= 0.0) smc(i,j,l) = spval + enddo + enddo + enddo + endif + + ! soill + if(trim(fieldname)=='soill') then + !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,sh2o,arrayr43d,sm,fillvalue,spval) + do l=1,nsoil + do j=jsta,jend + do i=ista, iend + sh2o(i,j,l) = arrayr43d(i,j,l) + if( abs(arrayr43d(i,j,l)-fillValue) < small) sh2o(i,j,l) = spval + if (sm(i,j) /= 0.0) sh2o(i,j,l) = spval + enddo + enddo + enddo + endif + ! model level ozone mixing ratio #ifdef MULTI_GASES if(trim(fieldname)=='spo3') then From 2d342231225b6d00875d164cb86e09ea335fe6b1 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 23 Feb 2024 11:54:56 -0700 Subject: [PATCH 5/7] Submodule pointer updates for ccpp-framework and ccpp-physics, fix compile error on macOS (#747) * Make ozone physics CCPP compliant by removing 'optional' and 'pointer' attributes: ccpp/physics PRs #168 #169 #150 * Bug fix in io/module_wrt_grid_comp.F90 to compile on macOS with apple-clang * Update ccpp-framework --- ccpp/framework | 2 +- ccpp/physics | 2 +- io/module_wrt_grid_comp.F90 | 2 ++ 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/ccpp/framework b/ccpp/framework index 219f2e9c8..f0b9a18b0 160000 --- a/ccpp/framework +++ b/ccpp/framework @@ -1 +1 @@ -Subproject commit 219f2e9c88b7b774becac2bd1453696e105af1c4 +Subproject commit f0b9a18b005d950cb9b0038fbc827b6b37500f43 diff --git a/ccpp/physics b/ccpp/physics index 51452b8b5..4ac3009e5 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 51452b8b52d5f3ded391847aefc21b086bb27b64 +Subproject commit 4ac3009e5181756a02acf95b5305ba016cd9a9cd diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 9a73f3d43..de0cedb6f 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -28,7 +28,9 @@ module module_wrt_grid_comp ! use mpi use esmf + use fms_mod, only : uppercase use fms + use mpp_mod, only : mpp_init, mpp_error use write_internal_state use module_fv3_io_def, only : num_pes_fcst, & From 4941026c5881ed623b39934905d2e4ecd9907b8a Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 29 Feb 2024 18:12:55 -0500 Subject: [PATCH 6/7] tisfc bugfix (#783) * ice model to determine the ice temperature over both sea ice and lake ice --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 4ac3009e5..cc114f40b 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 4ac3009e5181756a02acf95b5305ba016cd9a9cd +Subproject commit cc114f40bef4ca8d19ccd739cbaec3fc829a607c From caa61fc048a5eb8b72e8c988dea0883110f5e0ec Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 29 Mar 2024 12:41:14 -0400 Subject: [PATCH 7/7] update ccpp/physics commit hash --- ccpp/physics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ccpp/physics b/ccpp/physics index 1541f8f7b..b2a733836 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 1541f8f7b787f443fdd8b6812255f8492751ca00 +Subproject commit b2a733836e4a335a817120baec915fd21ea64f76