Skip to content

Stochastic physics changes #70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
MingjingTong-NOAA wants to merge 4 commits into
base: main
Choose a base branch
from
+39 −32
Open

Stochastic physics changes #70

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6be9fe4
Made changes to be able to use updated stochastic physics package (ne…
MingjingTong-NOAA Mar 5, 2026
7996c79
Revert a change
MingjingTong-NOAA Mar 5, 2026
b180252
Add iau_offset into atmos_model_nml and assign the value of Atmos%iau…
MingjingTong-NOAA Mar 12, 2026
b3c7720
Remove MPI communicator from Init_parm and passing the number of thre…
MingjingTong-NOAA Mar 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 39 additions & 32 deletions SHiELD/atmos_model.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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.
Copy link
Copy Markdown
Contributor

@JosephMouallem JosephMouallem Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a specific reason for this name change? I would prefer keeping the current name to avoid having to update all of our runscripts.

Copy link
Copy Markdown
Author

@MingjingTong-NOAA MingjingTong-NOAA Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I change the name, because I introduced another 'first_time_step' variable in physics here

./GFS_layer/GFS_typedefs.F90: logical :: first_time_step !< flag signaling first time step for time integration routine

This variable is used in the atmos_model.F90

./atmos_model.F90: IPD_control%first_time_step = .false.
./atmos_model.F90: IPD_control%first_time_step = .true.

Also, I think write_first_time_step better indicates the meaning of this variables.

logical :: fprint = .true.
logical :: ignore_rst_cksum = .false. ! enforce (.false.) or override (.true.) data integrity restart checksums

Expand All @@ -296,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, first_time_step, fdiag, fprint, &
fdiag_override, ignore_rst_cksum, fullcoupler_fluxes
blocksize, chksum_debug, dycore_only, debug, sync, write_first_time_step, fdiag, fprint, &
fdiag_override, ignore_rst_cksum, fullcoupler_fluxes, iau_offset

type (time_type) :: diag_time, diag_time_fhzero
logical :: fdiag_fix = .false.
Expand Down Expand Up @@ -519,6 +518,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
Expand Down Expand Up @@ -594,7 +596,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)
Expand Down Expand Up @@ -639,9 +641,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, nthrds, ierr)
if (ierr/=0) call mpp_error(FATAL, 'Call to stochastic_physics_wrapper failed')
endif
#endif

call mpp_clock_end(setupClock)
Expand Down Expand Up @@ -746,7 +750,6 @@ end subroutine atm_stock_pe
! </INOUT>

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
Expand All @@ -756,9 +759,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step, do_concurrent_ra
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
integer :: mlon, mlat, nlon, nlat, nlev, sec, dt, sec_prev
integer :: ierr, io, logunit
Expand Down Expand Up @@ -787,7 +788,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
Expand All @@ -805,6 +805,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)
Expand Down Expand Up @@ -883,29 +884,14 @@ 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
Init_parm%fn_nml='using internal file'

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
Expand Down Expand Up @@ -948,6 +934,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, nthrds, 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
Expand Down Expand Up @@ -1006,6 +1001,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)
Expand Down Expand Up @@ -1102,14 +1101,14 @@ 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(Atmos%iau_offset > zero) then
if( time_int - Atmos%iau_offset*3600. > zero ) then
time_int = time_int - Atmos%iau_offset*3600.
Expand Down Expand Up @@ -1373,6 +1372,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)
Expand Down