Skip to content

Commit

Permalink
Update gsl/develop from develop 2020/12/08 (NCAR#61)
Browse files Browse the repository at this point in the history 
* Fix for updating stochastic physics on separate time-step. (NCAR#199)

This bug fix allows the random patterns in the stochastic physics persist the for a period of time (defined as SKEBINT,SPPTINT, etc.) before calculating new patterns.
The fix is to move the allocation of the saved variables into the init section of stochastic_physics_wrapper, and remove the deallocates in the run section.

* Bug fixes in (1) running with frac_grid=T and GFDL MP and (2) restarting with frac_grid=T (NCAR#204)

* -- Pointing to Moorthi's modifications in ccpp/physics, which fixed the crash when running GFDL MP with frac_grid=T;
-- Not setting fice to zero in order to leave lake ice untouched;
-- Restart in the coupled model with the default physics is reproducible, if bad water temperature is only filtered at initial time;
Co-authored-with: Shrinivas Moorthi <[email protected]>
Co-authored-with: Denise Worthen <[email protected]>

* Revert change to .gitmodules and update submodule pointer for ccpp-physics

Co-authored-by: Phil Pegion <[email protected]>
Co-authored-by: shansun6 <[email protected]>
  • Loading branch information
@climbfuji @pjpegion @ShanSunNOAA
3 people authored Dec 8, 2020
1 parent bd71c2a commit c59787a
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 69 deletions.
8 changes: 5 additions & 3 deletions atmos_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ module atmos_model_mod
use IPD_driver, only: IPD_initialize, IPD_initialize_rst
use CCPP_driver, only: CCPP_step, non_uniform_blocks

use stochastic_physics_wrapper_mod, only: stochastic_physics_wrapper
use stochastic_physics_wrapper_mod, only: stochastic_physics_wrapper,stochastic_physics_wrapper_end
#else
use IPD_driver, only: IPD_initialize, IPD_initialize_rst, IPD_step
use physics_abstraction_layer, only: time_vary_step, radiation_step1, physics_step1, physics_step2
Expand Down Expand Up @@ -962,6 +962,9 @@ subroutine atmos_model_end (Atmos)
!---- termination routine for atmospheric model ----

call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst)

call stochastic_physics_wrapper_end(IPD_Control)

if(restart_endfcst) then
call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, &
IPD_Control, Atmos%domain)
Expand Down Expand Up @@ -1761,7 +1764,6 @@ subroutine assign_importdata(rc)
nb = Atm_block%blkno(i,j)
ix = Atm_block%ixp(i,j)

IPD_Data(nb)%Sfcprop%fice(ix) = zero
IPD_Data(nb)%Coupling%slimskin_cpl(ix) = IPD_Data(nb)%Sfcprop%slmsk(ix)
ofrac = IPD_Data(nb)%Sfcprop%oceanfrac(ix)
if (ofrac > zero) then
Expand All @@ -1776,7 +1778,7 @@ subroutine assign_importdata(rc)
if (abs(one-ofrac) < epsln) then
IPD_Data(nb)%Sfcprop%slmsk(ix) = zero
IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero
end if
endif
endif
endif
enddo
Expand Down
2 changes: 1 addition & 1 deletion ccpp/physics
Submodule physics updated 2 files
+94 −63 physics/GFS_surface_composites.F90
+27 −10 physics/GFS_surface_composites.meta
45 changes: 27 additions & 18 deletions io/FV3GFS_io.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1070,14 +1070,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
endif

if(Model%frac_grid) then ! obtain slmsk from landfrac
!! next 5 lines are temporary till lake model is available
if (Sfcprop(nb)%lakefrac(ix) > zero) then
! Sfcprop(nb)%lakefrac(ix) = nint(Sfcprop(nb)%lakefrac(ix))
Sfcprop(nb)%landfrac(ix) = one - Sfcprop(nb)%lakefrac(ix)
if (Sfcprop(nb)%lakefrac(ix) == zero) Sfcprop(nb)%fice(ix) = zero
endif
Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix))
if (Sfcprop(nb)%fice(ix) > Model%min_lakeice .and. Sfcprop(nb)%landfrac(ix) == zero) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates ice if co-exist
Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) !nint/floor are options
else ! obtain landfrac from slmsk
if (Sfcprop(nb)%slmsk(ix) > 1.9_r8) then
Sfcprop(nb)%landfrac(ix) = zero
Expand All @@ -1088,16 +1081,32 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)

if (Sfcprop(nb)%lakefrac(ix) > zero) then
Sfcprop(nb)%oceanfrac(ix) = zero ! lake & ocean don't coexist in a cell
! if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) then
! Sfcprop(nb)%fice(ix) = zero
! if (Sfcprop(nb)%slmsk(ix) == 2) Sfcprop(nb)%slmsk(ix) = 0
! endif
if (Sfcprop(nb)%slmsk(ix) /= one) then
if (Sfcprop(nb)%fice(ix) >= Model%min_lakeice) then
if (Sfcprop(nb)%slmsk(ix) < 1.9_r8) &
write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=2 at nb,ix=' &
,nb,ix,Sfcprop(nb)%fice(ix),Sfcprop(nb)%slmsk(ix),Sfcprop(nb)%lakefrac(ix)
Sfcprop(nb)%slmsk(ix) = 2.
else if (Sfcprop(nb)%slmsk(ix) > 1.e-7) then
write(*,'(a,2i3,3f6.2)') 'reset lake slmsk=0 at nb,ix=' &
,nb,ix,Sfcprop(nb)%fice(ix),Sfcprop(nb)%slmsk(ix),Sfcprop(nb)%lakefrac(ix)
Sfcprop(nb)%slmsk(ix) = zero
end if
end if
else
Sfcprop(nb)%oceanfrac(ix) = one - Sfcprop(nb)%landfrac(ix)
! if (Sfcprop(nb)%fice(ix) < Model%min_seaice) then
! Sfcprop(nb)%fice(ix) = zero
! if (Sfcprop(nb)%slmsk(ix) == 2) Sfcprop(nb)%slmsk(ix) = 0
! endif
if (Sfcprop(nb)%slmsk(ix) /= one) then
if (Sfcprop(nb)%fice(ix) >= Model%min_seaice) then
if (Sfcprop(nb)%slmsk(ix) < 1.9_r8) &
write(*,'(a,2i3,3f6.2)') 'reset sea slmsk=2 at nb,ix=' &
,nb,ix,Sfcprop(nb)%fice(ix),Sfcprop(nb)%slmsk(ix),Sfcprop(nb)%landfrac(ix)
Sfcprop(nb)%slmsk(ix) = 2.
else if (Sfcprop(nb)%slmsk(ix) > 1.e-7) then
write(*,'(a,2i3,4f6.2)') 'reset sea slmsk=0 at nb,ix=' &
,nb,ix,Sfcprop(nb)%fice(ix),Sfcprop(nb)%slmsk(ix),Sfcprop(nb)%landfrac(ix)
Sfcprop(nb)%slmsk(ix) = zero
end if
end if
endif
!
!--- NSSTM variables
Expand Down Expand Up @@ -1351,7 +1360,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
endif

if (sfc_var2(i,j,nvar_s2m) < -9990.0_r8) then
if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorli')
if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorlw')
!$omp parallel do default(shared) private(nb, ix)
do nb = 1, Atm_block%nblks
do ix = 1, Atm_block%blksz(nb)
Expand All @@ -1366,7 +1375,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
!$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))
if( Model%phour < 1.e-7) 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)%zorl(ix) = Sfcprop(nb)%zorll(ix) * Sfcprop(nb)%landfrac(ix) &
Expand Down
111 changes: 64 additions & 47 deletions stochastic_physics/stochastic_physics_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,25 +92,43 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
return
endif
end if
allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
if (GFS_Control%do_sppt) then
allocate(sppt_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if (GFS_Control%do_shum) then
allocate(shum_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if (GFS_Control%do_skeb) then
allocate(skebu_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
allocate(skebv_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast
allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%n_var_lndp))
end if
if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(stype(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(vfrac(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
endif

do nb=1,Atm_block%nblks
xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:)
xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:)
end do

if ( GFS_Control%lndp_type .EQ. 1 ) then ! this scheme sets perts once
! Copy blocked data into contiguous arrays; no need to copy sfc_wts in (intent out)
allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%n_var_lndp))
do nb=1,Atm_block%nblks
xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:)
xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:)
end do
call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, &
sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, &
nthreads=nthreads)
! Copy contiguous data back; no need to copy xlat/xlon, these are intent(in) - just deallocate
do nb=1,Atm_block%nblks
GFS_Data(nb)%Coupling%sfc_wts(:,:) = sfc_wts(nb,1:GFS_Control%blksz(nb),:)
end do
deallocate(xlat)
deallocate(xlon)
deallocate(sfc_wts)
end if
! Consistency check for cellular automata
Expand All @@ -126,27 +144,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
else initalize_stochastic_physics

if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .EQ. 2) ) then
! Copy blocked data into contiguous arrays; no need to copy weights in (intent(out))
allocate(xlat(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(xlon(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
do nb=1,Atm_block%nblks
xlat(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlat(:)
xlon(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Grid%xlon(:)
end do
if (GFS_Control%do_sppt) then
allocate(sppt_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if (GFS_Control%do_shum) then
allocate(shum_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if (GFS_Control%do_skeb) then
allocate(skebu_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
allocate(skebv_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%levs))
end if
if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast
allocate(sfc_wts(1:Atm_block%nblks,maxval(GFS_Control%blksz),1:GFS_Control%n_var_lndp))
end if

call run_stochastic_physics(GFS_Control%levs, GFS_Control%kdt, GFS_Control%phour, GFS_Control%blksz, xlat=xlat, xlon=xlon, &
sppt_wts=sppt_wts, shum_wts=shum_wts, skebu_wts=skebu_wts, skebv_wts=skebv_wts, sfc_wts=sfc_wts, &
nthreads=nthreads)
Expand All @@ -155,32 +152,23 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
do nb=1,Atm_block%nblks
GFS_Data(nb)%Coupling%sppt_wts(:,:) = sppt_wts(nb,1:GFS_Control%blksz(nb),:)
end do
deallocate(sppt_wts)
end if
if (GFS_Control%do_shum) then
do nb=1,Atm_block%nblks
GFS_Data(nb)%Coupling%shum_wts(:,:) = shum_wts(nb,1:GFS_Control%blksz(nb),:)
end do
deallocate(shum_wts)
end if
if (GFS_Control%do_skeb) then
do nb=1,Atm_block%nblks
GFS_Data(nb)%Coupling%skebu_wts(:,:) = skebu_wts(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Coupling%skebv_wts(:,:) = skebv_wts(nb,1:GFS_Control%blksz(nb),:)
end do
deallocate(skebu_wts)
deallocate(skebv_wts)
end if
if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
do nb=1,Atm_block%nblks
GFS_Data(nb)%Coupling%sfc_wts(:,:) = sfc_wts(nb,1:GFS_Control%blksz(nb),:)
end do

allocate(smc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(slc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(stc(1:Atm_block%nblks,maxval(GFS_Control%blksz),GFS_Control%lsoil))
allocate(stype(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
allocate(vfrac(1:Atm_block%nblks,maxval(GFS_Control%blksz)))
do nb=1,Atm_block%nblks
stype(nb,1:GFS_Control%blksz(nb)) = GFS_Data(nb)%Sfcprop%stype(:)
smc(nb,1:GFS_Control%blksz(nb),:) = GFS_Data(nb)%Sfcprop%smc(:,:)
Expand All @@ -202,21 +190,13 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
write(6,*) 'call to GFS_apply_lndp failed'
return
endif
deallocate(stype)
deallocate(sfc_wts)
do nb=1,Atm_block%nblks
GFS_Data(nb)%Sfcprop%smc(:,:) = smc(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Sfcprop%slc(:,:) = slc(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Sfcprop%stc(:,:) = stc(nb,1:GFS_Control%blksz(nb),:)
GFS_Data(nb)%Sfcprop%vfrac(:) = vfrac(nb,1:GFS_Control%blksz(nb))
enddo
deallocate(smc)
deallocate(slc)
deallocate(stc)
deallocate(vfrac)
endif ! lndp block
deallocate(xlat)
deallocate(xlon)
end if

endif initalize_stochastic_physics
Expand Down Expand Up @@ -309,4 +289,41 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)

end subroutine stochastic_physics_wrapper


subroutine stochastic_physics_wrapper_end (GFS_Control)

use GFS_typedefs, only: GFS_control_type, GFS_data_type
use stochastic_physics, only: finalize_stochastic_physics

implicit none

type(GFS_control_type), intent(inout) :: GFS_Control

if (GFS_Control%do_sppt .OR. GFS_Control%do_shum .OR. GFS_Control%do_skeb .OR. (GFS_Control%lndp_type .GT. 0) ) then
if (allocated(xlat)) deallocate(xlat)
if (allocated(xlon)) deallocate(xlon)
if (GFS_Control%do_sppt) then
if (allocated(sppt_wts)) deallocate(sppt_wts)
end if
if (GFS_Control%do_shum) then
if (allocated(shum_wts)) deallocate(shum_wts)
end if
if (GFS_Control%do_skeb) then
if (allocated(skebu_wts)) deallocate(skebu_wts)
if (allocated(skebv_wts)) deallocate(skebv_wts)
end if
if ( GFS_Control%lndp_type .EQ. 2 ) then ! this scheme updates through forecast
if (allocated(sfc_wts)) deallocate(sfc_wts)
end if
if (GFS_Control%lndp_type .EQ. 2) then ! save wts, and apply lndp scheme
if (allocated(smc)) deallocate(smc)
if (allocated(slc)) deallocate(slc)
if (allocated(stc)) deallocate(stc)
if (allocated(stype)) deallocate(stype)
if (allocated(vfrac)) deallocate(vfrac)
endif
call finalize_stochastic_physics()
endif
end subroutine stochastic_physics_wrapper_end

end module stochastic_physics_wrapper_mod

0 comments on commit c59787a

Please sign in to comment.