From 6be9fe49c71d338e308c8fc52f22b4fa207292b0 Mon Sep 17 00:00:00 2001 From: "Mingjing.Tong" Date: Thu, 5 Mar 2026 17:24:28 -0500 Subject: [PATCH 1/4] Made changes to be able to use updated stochastic physics package (new repo) --- SHiELD/atmos_model.F90 | 70 ++++++++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 30 deletions(-) diff --git a/SHiELD/atmos_model.F90 b/SHiELD/atmos_model.F90 index 9e65de3..4c3a5aa 100644 --- a/SHiELD/atmos_model.F90 +++ b/SHiELD/atmos_model.F90 @@ -89,9 +89,7 @@ module atmos_model_mod use diag_integral_mod, only : diag_integral_init #ifdef STOCHY -use stochastic_physics, only: init_stochastic_physics, & - run_stochastic_physics -use stochastic_physics_sfc, only: run_stochastic_physics_sfc +use stochastic_physics_wrapper_mod, only: stochastic_physics_wrapper,stochastic_physics_wrapper_end #endif use FV3GFS_io_mod, only: FV3GFS_restart_read, FV3GFS_restart_write, & FV3GFS_IPD_checksum, & @@ -280,7 +278,7 @@ module atmos_model_mod logical :: dycore_only = .false. logical :: debug = .false. logical :: sync = .false. -logical :: first_time_step = .false. +logical :: write_first_time_step = .false. logical :: fprint = .true. logical :: ignore_rst_cksum = .false. ! enforce (.false.) or override (.true.) data integrity restart checksums @@ -297,7 +295,7 @@ module atmos_model_mod integer :: nxblocks = 1 integer :: nyblocks = 1 namelist /atmos_model_nml/ do_netcdf_restart, restart_tbot_qbot, nxblocks, nyblocks, & - blocksize, chksum_debug, dycore_only, debug, sync, first_time_step, fdiag, fprint, & + blocksize, chksum_debug, dycore_only, debug, sync, write_first_time_step, fdiag, fprint, & fdiag_override, ignore_rst_cksum, fullcoupler_fluxes type (time_type) :: diag_time, diag_time_fhzero @@ -519,6 +517,9 @@ subroutine update_atmos_model_up( Surface_boundary, Atmos ) call getiauforcing(IPD_Control,IAU_data) if (mpp_pe() == mpp_root_pe() .and. debug) write(6,*) "end of radiation and physics step" + ! Update flag for first time step of time integration + IPD_control%first_time_step = .false. + call mpp_clock_end(shieldClock) ! Reupdate the phys fluxes here, as the coupler fluxes over land points get reupdated @@ -594,7 +595,7 @@ subroutine update_atmos_model_radiation (Surface_boundary, Atmos) ! name change ! AM4/src/atmos_phys/atmos_param/radiation/driver/radiation_driver.F90 !--- local variables--- integer :: nb, jdat(8) - integer :: nthrds + integer :: nthrds, ierr call set_atmosphere_pelist() ! should be called before local clocks since they are defined on local atm(n)%pelist call mpp_clock_begin(shieldClock) @@ -639,9 +640,11 @@ subroutine update_atmos_model_radiation (Surface_boundary, Atmos) ! name change #ifdef STOCHY !--- call stochastic physics pattern generation / cellular automata - if (IPD_Control%do_sppt .OR. IPD_Control%do_shum .OR. IPD_Control%do_skeb .OR. IPD_Control%do_sfcperts) then - call run_stochastic_physics(IPD_Control, IPD_Data(:)%Grid, IPD_Data(:)%Coupling, nthrds) - end if + if (IPD_Control%do_sppt .or. IPD_Control%do_shum .or. IPD_Control%do_skeb .or. & + IPD_Control%lndp_type > 0 .or. IPD_Control%do_ca .or. IPD_Control%do_spp) then + call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, ierr) + if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') + endif #endif call mpp_clock_end(setupClock) @@ -746,17 +749,17 @@ end subroutine atm_stock_pe ! subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_radiation) !argument change to match the full coupler -!subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, iau_offset) !check how iau_offset should work #ifdef OPENMP use omp_lib #endif use mpp_mod, only: mpp_npes +#include "mpif.h" + type (atmos_data_type), intent(inout) :: Atmos type (time_type), intent(in) :: Time_init, Time, Time_step logical, intent(in) :: do_concurrent_radiation - !integer, intent(inout) :: iau_offset = 0 !--- local variables --- integer :: iau_offset = 0 integer :: unit, ntdiag, ntfamily, i, j, k @@ -787,7 +790,6 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra Atmos % Time_init = Time_init Atmos % Time = Time Atmos % Time_step = Time_step - Atmos % iau_offset = iau_offset call get_time (Atmos % Time_step, sec) call get_time (Atmos%Time - Atmos%Time_init, sec_prev) dt_phys = real(sec) ! integer seconds @@ -856,6 +858,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra !--- setup IPD Init_parm Init_parm%me = mpp_pe() Init_parm%master = mpp_root_pe() + Init_parm%fcst_mpi_comm%mpi_val = MPI_COMM_WORLD Init_parm%tile_num = tile_num Init_parm%isc = isc Init_parm%jsc = jsc @@ -883,6 +886,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra Init_parm%hydro = hydro Init_parm%do_inline_mp = do_inline_mp Init_parm%do_cosp = do_cosp + Init_parm%restart = Atm(mygrid)%flagstruct%warm_start allocate(Init_parm%input_nml_file, mold=input_nml_file) Init_parm%input_nml_file => input_nml_file @@ -890,22 +894,6 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra call IPD_initialize (IPD_Control, IPD_Data, IPD_Diag, IPD_Restart, Init_parm) - IPD_Control%sfc_coupled = fullcoupler_fluxes - -#ifdef STOCHY - if (IPD_Control%do_sppt .OR. IPD_Control%do_shum .OR. IPD_Control%do_skeb .OR. IPD_Control%do_sfcperts) then - ! Initialize stochastic physics - call init_stochastic_physics(IPD_Control, Init_parm, mpp_npes(), nthrds) - if (mpp_pe() == mpp_root_pe()) print *,'do_skeb=',IPD_Control%do_skeb - end if - - if (IPD_Control%do_sfcperts) then - ! Get land surface perturbations here (move to GFS_time_vary - ! step if wanting to update each time-step) - call run_stochastic_physics_sfc(IPD_Control, IPD_Data(:)%Grid, IPD_Data(:)%Coupling) - end if -#endif - Atm(mygrid)%flagstruct%do_diss_est = IPD_Control%do_skeb ! initialize the IAU module @@ -948,6 +936,15 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra call FV3GFS_IPD_checksum(IPD_Control, IPD_Data, Atm_block) endif +#ifdef STOCHY + if (IPD_Control%do_sppt .or. IPD_Control%do_shum .or. IPD_Control%do_skeb .or. & + IPD_Control%lndp_type > 0 .or. IPD_Control%do_ca .or. IPD_Control%do_spp) then +!--- Initialize stochastic physics pattern generation / cellular automata for first time step + call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, ierr) + if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') + end if +#endif + !--- set the initial diagnostic timestamp diag_time = Time if (Atmos%iau_offset > zero) then @@ -1006,6 +1003,10 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra endif shieldClock= mpp_clock_id( '--SHiELD Physics ', flags=clock_flag_default, grain=CLOCK_COMPONENT ) + ! Set flag for first time step of time integration + IPD_control%first_time_step = .true. + +!----------------------------------------------------------------------- call diag_integral_init (Atmos % Time_init, Atmos % Time, & Atmos % lon_bnd(:,:), & Atmos % lat_bnd(:,:), Atmos % area) @@ -1102,14 +1103,15 @@ subroutine update_atmos_model_state (Atmos) time_int = real(isec) ! For overflow safety cast scaled fdiag values to 64-bit nearest integers ! before interacting with seconds. - if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. first_time_step) ) then + if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. write_first_time_step) ) then if (mpp_pe() == mpp_root_pe()) write(6,*) "---isec,seconds",isec,seconds if (mpp_pe() == mpp_root_pe()) write(6,*) ' gfs diags time since last bucket empty: ',time_int/3600.,'hrs' call atmosphere_nggps_diag(Atmos%Time) endif ! For overflow safety cast scaled fdiag values to 64-bit nearest integers ! before interacting with seconds. - if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. first_time_step)) then + !if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. write_first_time_step)) then + if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. write_first_time_step) then if(Atmos%iau_offset > zero) then if( time_int - Atmos%iau_offset*3600. > zero ) then time_int = time_int - Atmos%iau_offset*3600. @@ -1373,6 +1375,14 @@ subroutine atmos_model_end (Atmos) !---- termination routine for atmospheric model ---- call atmosphere_end (Atmos % Time, Atmos%grid) + +#ifdef STOCHY + if (IPD_Control%do_sppt .or. IPD_Control%do_shum .or. IPD_Control%do_skeb .or. & + IPD_Control%lndp_type > 0 .or. IPD_Control%do_ca .or. IPD_Control%do_spp) then + call stochastic_physics_wrapper_end(IPD_control) + endif +#endif + if (.not. dycore_only) then call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, & IPD_Control, Atmos%domain) From 7996c797b16f2b52f2faf53f68e623d8966773f9 Mon Sep 17 00:00:00 2001 From: "Mingjing.Tong" Date: Thu, 5 Mar 2026 17:31:29 -0500 Subject: [PATCH 2/4] Revert a change --- SHiELD/atmos_model.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SHiELD/atmos_model.F90 b/SHiELD/atmos_model.F90 index 4c3a5aa..392db63 100644 --- a/SHiELD/atmos_model.F90 +++ b/SHiELD/atmos_model.F90 @@ -1110,8 +1110,7 @@ subroutine update_atmos_model_state (Atmos) endif ! For overflow safety cast scaled fdiag values to 64-bit nearest integers ! before interacting with seconds. - !if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. write_first_time_step)) then - if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. write_first_time_step) then + if (ANY(nint(fdiag(:)*3600.0, kind=8) == seconds) .or. (fdiag_fix .and. mod(seconds, nint(fdiag(1)*3600.0, kind=8)) .eq. 0) .or. (IPD_Control%kdt == 1 .and. write_first_time_step)) then if(Atmos%iau_offset > zero) then if( time_int - Atmos%iau_offset*3600. > zero ) then time_int = time_int - Atmos%iau_offset*3600. From b180252ca5c33e6294465e3e8cff8c976c76ad40 Mon Sep 17 00:00:00 2001 From: "Mingjing.Tong" Date: Thu, 12 Mar 2026 17:54:28 -0400 Subject: [PATCH 3/4] Add iau_offset into atmos_model_nml and assign the value of Atmos%iau_offset in SHiELD/atmos_model.F90 --- SHiELD/atmos_model.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SHiELD/atmos_model.F90 b/SHiELD/atmos_model.F90 index 392db63..57caa3d 100644 --- a/SHiELD/atmos_model.F90 +++ b/SHiELD/atmos_model.F90 @@ -294,9 +294,10 @@ module atmos_model_mod logical :: restart_tbot_qbot = .false. integer :: nxblocks = 1 integer :: nyblocks = 1 +integer :: iau_offset = 0 namelist /atmos_model_nml/ do_netcdf_restart, restart_tbot_qbot, nxblocks, nyblocks, & blocksize, chksum_debug, dycore_only, debug, sync, write_first_time_step, fdiag, fprint, & - fdiag_override, ignore_rst_cksum, fullcoupler_fluxes + fdiag_override, ignore_rst_cksum, fullcoupler_fluxes, iau_offset type (time_type) :: diag_time, diag_time_fhzero logical :: fdiag_fix = .false. @@ -761,7 +762,6 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra type (time_type), intent(in) :: Time_init, Time, Time_step logical, intent(in) :: do_concurrent_radiation !--- local variables --- - integer :: iau_offset = 0 integer :: unit, ntdiag, ntfamily, i, j, k integer :: mlon, mlat, nlon, nlat, nlev, sec, dt, sec_prev integer :: ierr, io, logunit @@ -807,6 +807,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra IF ( file_exists('input.nml')) THEN read(input_nml_file, nml=atmos_model_nml, iostat=io) ierr = fms_check_nml_error(io, 'atmos_model_nml') + Atmos%iau_offset = iau_offset endif !----------------------------------------------------------------------- call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) From b3c772047a1c36d07a9bdd48d0f01147fea330fb Mon Sep 17 00:00:00 2001 From: Mingjing Tong Date: Sun, 29 Mar 2026 14:19:09 -0400 Subject: [PATCH 4/4] Remove MPI communicator from Init_parm and passing the number of threads into stochastic_physics_wrapper --- SHiELD/atmos_model.F90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/SHiELD/atmos_model.F90 b/SHiELD/atmos_model.F90 index 57caa3d..1fbe41b 100644 --- a/SHiELD/atmos_model.F90 +++ b/SHiELD/atmos_model.F90 @@ -643,7 +643,7 @@ subroutine update_atmos_model_radiation (Surface_boundary, Atmos) ! name change !--- call stochastic physics pattern generation / cellular automata if (IPD_Control%do_sppt .or. IPD_Control%do_shum .or. IPD_Control%do_skeb .or. & IPD_Control%lndp_type > 0 .or. IPD_Control%do_ca .or. IPD_Control%do_spp) then - call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, ierr) + call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, nthrds, ierr) if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') endif #endif @@ -756,8 +756,6 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra #endif use mpp_mod, only: mpp_npes -#include "mpif.h" - type (atmos_data_type), intent(inout) :: Atmos type (time_type), intent(in) :: Time_init, Time, Time_step logical, intent(in) :: do_concurrent_radiation @@ -859,7 +857,6 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra !--- setup IPD Init_parm Init_parm%me = mpp_pe() Init_parm%master = mpp_root_pe() - Init_parm%fcst_mpi_comm%mpi_val = MPI_COMM_WORLD Init_parm%tile_num = tile_num Init_parm%isc = isc Init_parm%jsc = jsc @@ -941,7 +938,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra if (IPD_Control%do_sppt .or. IPD_Control%do_shum .or. IPD_Control%do_skeb .or. & IPD_Control%lndp_type > 0 .or. IPD_Control%do_ca .or. IPD_Control%do_spp) then !--- Initialize stochastic physics pattern generation / cellular automata for first time step - call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, ierr) + call stochastic_physics_wrapper(IPD_control, IPD_data, Atm_block, nthrds, ierr) if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed') end if #endif