From 25ed5ef652fe30284192c9eb6ae4ef5882516863 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Tue, 22 Dec 2020 18:40:13 +0000 Subject: [PATCH] additions for stochy restarts --- config_src/nuopc_driver/mom_cap.F90 | 11 ++++++++--- config_src/solo_driver/MOM_driver.F90 | 3 +++ src/core/MOM.F90 | 14 +++++++------- .../vertical/MOM_diabatic_driver.F90 | 16 +++++++++++++--- .../vertical/MOM_energetic_PBL.F90 | 19 ++++++++++--------- 5 files changed, 41 insertions(+), 22 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 6880addfb2..12d70fe4e7 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -91,7 +91,9 @@ module MOM_cap_mod use NUOPC_Model, only: model_label_SetRunClock => label_SetRunClock use NUOPC_Model, only: model_label_Finalize => label_Finalize use NUOPC_Model, only: SetVM +#ifdef UFS use get_stochy_pattern_mod, only: write_stoch_restart_ocn +#endif implicit none; private @@ -1587,14 +1589,17 @@ subroutine ModelAdvance(gcomp, rc) call ocean_model_restart(ocean_state, restartname=restartname) ! write stochastic physics restart file if active +#ifdef UFS if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then - write(restartname,'(A)')"ocn_stoch.res.nc") + write(restartname,'(A)')"ocn_stoch.res.nc" else write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') & - "oc_stoch.res.", year, month, day, hour, minute, seconds,".nc" + "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc" endif call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) - call write_stoch_restart_ocn('RESTART/'//trim(timestamp)//'.ocn_stoch.res.nc') + if (is_root_pe()) print*,'calling write_stoch_restart_ocn ',trim(restartname) + call write_stoch_restart_ocn('RESTART/'//trim(restartname)) +#endif endif if (is_root_pe()) then diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index ba52d9c02a..81fb79207f 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -74,6 +74,9 @@ program MOM_main use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves +#ifdef UFS + use get_stochy_pattern_mod, only: write_stoch_restart_ocn +#endif implicit none diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a4f2f81af2..5143ffbf77 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -142,7 +142,9 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline +#ifdef UFS use stochastic_physics, only : init_stochastic_physics_ocn +#endif implicit none ; private @@ -229,8 +231,6 @@ module MOM logical :: offline_tracer_mode = .false. !< If true, step_offline() is called instead of step_MOM(). !! This is intended for running MOM6 in offline tracer mode - logical :: do_stochy = .false. - !< If true, call stochastic physics pattern generator type(time_type), pointer :: Time !< pointer to the ocean clock real :: dt !< (baroclinic) dynamics time step [T ~> s] @@ -2361,6 +2361,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) call destroy_dyn_horgrid(dG_in) + do_epbl=.false. + do_sppt=.false. +#ifdef UFS num_procs=num_PEs() allocate(pelist(num_procs)) call Get_PElist(pelist,commID = mom_comm) @@ -2368,12 +2371,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & master=root_PE() !print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) - do_epbl=.false. - do_sppt=.false. + if (master) print*,'about to call init_stochastic_physics' call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,do_epbl,do_sppt,master,mom_comm,iret) - if (do_sppt .eq. .true.) CS%do_stochy=.true. - if (do_epbl .eq. .true.) CS%do_stochy=.true. - !print*,'back from init_stochastic_physics_ocn',CS%do_stochy +#endif ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 15d41293c1..78068ba44d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -69,7 +69,9 @@ module MOM_diabatic_driver use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_wave_speed, only : wave_speeds use MOM_wave_interface, only : wave_parameters_CS +#ifdef UFS use stochastic_physics, only : run_stochastic_physics_ocn +#endif implicit none ; private @@ -292,23 +294,28 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & real, dimension(SZI_(G),SZJ_(G)) :: sppt_wts real, dimension(SZI_(G),SZJ_(G),2) :: t_rp +#ifdef UFS real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_in real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_in !< thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_in !< thickness [H ~> m or kg m-2] real :: t_tend,s_tend,h_tend ! holder for tendencey needed for SPPT real :: t_pert,s_pert,h_pert ! holder for tendencey needed for SPPT +#endif if (G%ke == 1) return +#ifdef UFS ! save copy of the date for SPPT if (CS%do_sppt) then h_in=h t_in=tv%T s_in=tv%S endif + print*,'calling run_stochastic_physics' call run_stochastic_physics_ocn(t_rp,sppt_wts) !print*,'in diabatic',CS%do_sppt,size(t_in,1),size(t_in,2),size(t_in,3),size(sppt_wts,1),size(sppt_wts,2) - !print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts) + print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts) + print*,'in diabatic',CS%do_sppt,minval(t_rp),maxval(t_rp) if (CS%id_t_rp1 > 0) then call post_data(CS%id_t_rp1, t_rp(:,:,1), CS%diag) endif @@ -318,6 +325,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_sppt_wts > 0) then call post_data(CS%id_sppt_wts, sppt_wts, CS%diag) endif +#endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -465,6 +473,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US) +#ifdef UFS if (CS%do_sppt) then do k=1,nz do j=js,je @@ -488,6 +497,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo enddo endif +#endif end subroutine diabatic @@ -891,7 +901,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, d endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -1623,7 +1633,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 277f714f2a..3c82049a8b 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -246,7 +246,7 @@ module MOM_energetic_PBL !! mixed layer model. It assumes that heating, cooling and freshwater fluxes !! have already been applied. All calculations are done implicitly, and there !! is no stability limit on the time step. -subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_int, G, GV, US, CS, & dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, & dT_expected, dS_expected, Waves ) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. @@ -285,6 +285,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, !! call to mixedlayer_init. real, dimension(SZI_(G),SZJ_(G),2), & intent(in) :: t_rp !< random pattern to perturb wind + logical, intent(in) :: stoch_epbl !< flag to pertrub production and dissipation of TKE real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]. real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less @@ -438,8 +439,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, do K=1,nz+1 ; Kd(K) = 0.0 ; enddo ! Make local copies of surface forcing and process them. - !print*,'PJP EPBL',minval(t_rp),maxval(t_rp) - u_star = fluxes%ustar(i,j)!*t_rp(i,j) + u_star = fluxes%ustar(i,j) u_star_Mean = fluxes%ustar_gustless(i,j) B_flux = buoy_flux(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then @@ -463,7 +463,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & - US, CS, eCD, t_rp(i,j,1),t_rp(i,j,2), dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) + US, CS, eCD, t_rp(i,j,1),t_rp(i,j,2), stoch_epbl, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) ! applly stochastic perturbation to TKE generation @@ -548,7 +548,7 @@ end subroutine energetic_PBL !! mixed layer model for a single column of water. subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, & - t_rp1,t_rp2, dt_diag, Waves, G, i, j) + t_rp1,t_rp2, stoch_epbl, dt_diag, Waves, G, i, j) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. @@ -587,7 +587,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !! call to mixedlayer_init. type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. real, intent(in) :: t_rp1 !< random value to perturb TKE production - real, intent(in) :: t_rp2 !< random value to perturb TKE production + real, intent(in) :: t_rp2 !< random value to perturb TKE dissipation + logical, intent(in) :: stoch_epbl !< flag to pertrub production and dissipation of TKE real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less !! than dt if there are two calls to mixedlayer [T ~> s]. type(wave_parameters_CS), & @@ -888,8 +889,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs else mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) endif - ! stochastically pertrub mech_TKE - mech_TKE=mech_TKE*t_rp1 + ! stochastically pertrub mech_TKE in the UFS + if (stoch_epbl) mech_TKE=mech_TKE*t_rp1 if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -973,7 +974,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs if (CS%TKE_diagnostics) & !eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag eCD%dTKE_mech_decay = exp_kh - mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) + if (stoch_epbl) mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) ! Accumulate any convectively released potential energy to contribute ! to wstar and to drive penetrating convection.