diff --git a/atmos_model.F90 b/atmos_model.F90 index 860079949..051f5918d 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/ccpp/physics b/ccpp/physics index bafcc9ebb..91b772665 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit bafcc9ebb6788696666e7592cf7577db81e5fbbe +Subproject commit 91b77266541d8a80dd23d09f80ee1e72e34af2d9 diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index a710941e2..46fdd8779 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -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 @@ -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 @@ -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) @@ -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) & diff --git a/stochastic_physics/stochastic_physics_wrapper.F90 b/stochastic_physics/stochastic_physics_wrapper.F90 index b5a1be065..5a3701aa8 100644 --- a/stochastic_physics/stochastic_physics_wrapper.F90 +++ b/stochastic_physics/stochastic_physics_wrapper.F90 @@ -92,16 +92,36 @@ 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) @@ -109,8 +129,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr) 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 @@ -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) @@ -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(:,:) @@ -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 @@ -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