Merge pull request #5 from MJHarrison-GFDL/ice_shelf_dynamic_sfc_accum

Ice shelf dynamic sfc accum

README.md

@@ -1,5 +1,5 @@
 [![Read The Docs Status](https://readthedocs.org/projects/mom6/badge/?badge=latest)](http://mom6.readthedocs.io/)
-[![codecov](https://codecov.io/gh/NOAA-GFDL/MOM6/branch/dev%2Fmaster/graph/badge.svg)](https://codecov.io/gh/NOAA-GFDL/MOM6)
+[![codecov](https://codecov.io/gh/NOAA-GFDL/MOM6/branch/dev/gfdl/graph/badge.svg?token=uF8SVydCdp)](https://codecov.io/gh/NOAA-GFDL/MOM6)
 
 # MOM6
 

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -113,6 +113,8 @@ module MOM_surface_forcing_gfdl
   real    :: Flux_const                     !< Piston velocity for surface restoring [Z T-1 ~> m s-1]
   real    :: Flux_const_salt                !< Piston velocity for surface salt restoring [Z T-1 ~> m s-1]
   real    :: Flux_const_temp                !< Piston velocity for surface temp restoring [Z T-1 ~> m s-1]
+  logical :: trestore_SPEAR_ECDA            !< If true, modify restoring data wrt local SSS
+  real    :: SPEAR_dTf_dS                   !< The derivative of the freezing temperature with salinity.
   logical :: salt_restore_as_sflux          !< If true, SSS restore as salt flux instead of water flux
   logical :: adjust_net_srestore_to_zero    !< Adjust srestore to zero (for both salt_flux or vprec)
   logical :: adjust_net_srestore_by_scaling !< Adjust srestore w/o moving zero contour
@@ -346,7 +348,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     open_ocn_mask(:,:) = 1.0
     if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice
       do j=js,je ; do i=is,ie
-        if (sfc_state%SST(i,j) <= -0.0539*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
+        if (sfc_state%SST(i,j) <= CS%SPEAR_dTf_dS*sfc_state%SSS(i,j)) open_ocn_mask(i,j)=0.0
       enddo ; enddo
     endif
     if (CS%salt_restore_as_sflux) then
@@ -400,6 +402,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   ! SST restoring logic
   if (CS%restore_temp) then
     call time_interp_external(CS%id_trestore, Time, data_restore)
+    if ( CS%trestore_SPEAR_ECDA ) then
+      do j=js,je ; do i=is,ie
+        if (abs(data_restore(i,j)+1.8)<0.0001) then
+          data_restore(i,j) = CS%SPEAR_dTf_dS*sfc_state%SSS(i,j)
+        endif
+      enddo ; enddo
+    endif
+
     do j=js,je ; do i=is,ie
       delta_sst = data_restore(i,j)- sfc_state%SST(i,j)
       delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore)
@@ -1448,7 +1458,15 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
                  "If true, read a file (temp_restore_mask) containing "//&
                  "a mask for SST restoring.", default=.false.)
 
+    call get_param(param_file, mdl, "SPEAR_ECDA_SST_RESTORE_TFREEZE", CS%trestore_SPEAR_ECDA, &
+                 "If true, modify SST restoring field using SSS state. This only modifies the "//&
+                 "restoring data that is within 0.0001degC of -1.8degC.", default=.false.)
+  else
+    CS%trestore_SPEAR_ECDA = .false. ! Needed to toggle logging of SPEAR_DTFREEZE_DS
   endif
+  call get_param(param_file, mdl, "SPEAR_DTFREEZE_DS", CS%SPEAR_dTf_dS, &
+                 "The derivative of the freezing temperature with salinity.", &
+                 units="deg C PSU-1", default=-0.054, do_not_log=.not.CS%trestore_SPEAR_ECDA)
 
   ! Optionally read tidal amplitude from input file [Z T-1 ~> m s-1] on model grid.
   ! Otherwise use default tidal amplitude for bottom frictionally-generated

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -97,6 +97,7 @@ 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
+
 !$use omp_lib             , only : omp_set_num_threads
 
 implicit none; private
@@ -1524,7 +1525,7 @@ subroutine ModelAdvance(gcomp, rc)
   integer                                :: nc
   type(ESMF_Time)                        :: MyTime
   integer                                :: seconds, day, year, month, hour, minute
-  character(ESMF_MAXSTR)                 :: restartname, cvalue
+  character(ESMF_MAXSTR)                 :: restartname, cvalue, stoch_restartname
   character(240)                         :: msgString
   character(ESMF_MAXSTR)                 :: casename
   integer                                :: iostat
@@ -1738,14 +1739,19 @@ subroutine ModelAdvance(gcomp, rc)
         ! write the final restart without a timestamp
         if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
            write(restartname,'(A)')"MOM.res"
+           write(stoch_restartname,'(A)')"ocn_stoch.res.nc"
         else
            write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2)') &
                 "MOM.res.", year, month, day, hour, minute, seconds
+           write(stoch_restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') &
+                "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc"
         endif
         call ESMF_LogWrite("MOM_cap: Writing restart :  "//trim(restartname), ESMF_LOGMSG_INFO)
 
         ! write restart file(s)
-        call ocean_model_restart(ocean_state, restartname=restartname)
+        call ocean_model_restart(ocean_state, restartname=restartname, &
+                                stoch_restartname=stoch_restartname)
+
      endif
 
      if (is_root_pe()) then

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -62,6 +62,7 @@ module MOM_ocean_model_nuopc
 use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum
 use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS
 use MOM_surface_forcing_nuopc, only : forcing_save_restart
+use get_stochy_pattern_mod,  only : write_stoch_restart_ocn
 use iso_fortran_env,           only : int64
 
 #include <MOM_memory.h>
@@ -176,6 +177,10 @@ module MOM_ocean_model_nuopc
                               !! steps can span multiple coupled time steps.
   logical :: diabatic_first   !< If true, apply diabatic and thermodynamic
                               !! processes before time stepping the dynamics.
+  logical :: do_sppt         !< If true, stochastically perturb the diabatic and
+                             !! write restarts
+  logical :: pert_epbl       !< If true, then randomly perturb the KE dissipation and
+                             !! genration termsand write restarts
 
   real :: eps_omesh           !< Max allowable difference between ESMF mesh and MOM6
                               !! domain coordinates
@@ -425,6 +430,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   endif
 
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
+! get number of processors and PE list for stocasthci physics initialization
+  call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, &
+                 "If true, then stochastically perturb the thermodynamic "//&
+                 "tendencies of T,S, and h.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  call get_param(param_file, mdl, "PERT_EPBL", OS%pert_epbl, &
+                 "If true, then stochastically perturb the kinetic energy "//&
+                 "production and dissipation terms.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
 
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
@@ -686,14 +702,17 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 end subroutine update_ocean_model
 
 !> This subroutine writes out the ocean model restart file.
-subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
+subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, num_rest_files)
   type(ocean_state_type),     pointer    :: OS !< A pointer to the structure containing the
                                                !! internal ocean state being saved to a restart file
   character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be
                                                !! prepended to the file name. (Currently this is unused.)
   character(len=*), optional, intent(in) :: restartname !< Name of restart file to use
                                                !! This option distinguishes the cesm interface from the
                                                !! non-cesm interface
+  character(len=*), optional, intent(in) :: stoch_restartname !< Name of restart file to use
+                                               !! This option distinguishes the cesm interface from the
+                                               !! non-cesm interface
   integer, optional, intent(out)         :: num_rest_files !< number of restart files written
 
   if (.not.MOM_state_is_synchronized(OS%MOM_CSp)) &
@@ -733,6 +752,11 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files)
         endif
      endif
   endif
+  if (present(stoch_restartname)) then
+      if (OS%do_sppt .OR. OS%pert_epbl) then
+         call write_stoch_restart_ocn('RESTART/'//trim(stoch_restartname))
+     endif
+  endif
 
 end subroutine ocean_model_restart
 ! </SUBROUTINE> NAME="ocean_model_restart"

config_src/drivers/solo_driver/MOM_driver.F90

@@ -332,7 +332,8 @@ program MOM_main
                  "The default value is given by DT.", units="s", default=dt)
   if (offline_tracer_mode) then
     call get_param(param_file, mod_name, "DT_OFFLINE", dt_forcing, &
-                   "Time step for the offline time step")
+                   "Length of time between reading in of input fields", &
+                   units='s', fail_if_missing=.true.)
     dt = dt_forcing
   endif
   ntstep = MAX(1,ceiling(dt_forcing/dt - 0.001))

config_src/external/GFDL_ocean_BGC/generic_tracer.F90

@@ -42,7 +42,7 @@ subroutine generic_tracer_init(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,axes,grid
     integer,                       intent(in) :: jsd !< Data start index in j direction
     integer,                       intent(in) :: jed !< Data end index in j direction
     integer,                       intent(in) :: nk  !< Number of levels in k direction
-    integer,                       intent(in) :: ntau !< Unknown
+    integer,                       intent(in) :: ntau !< The number of tracer time levels (always 1 for MOM6)
     integer,                       intent(in) :: axes(3) !< Domain axes?
     type(time_type),               intent(in) :: init_time !< Time
     real, dimension(:,:,:),target, intent(in) :: grid_tmask !< Mask
@@ -61,52 +61,61 @@ end subroutine generic_tracer_coupler_get
   !> Unknown
   subroutine  generic_tracer_coupler_accumulate(IOB_struc, weight, model_time)
     type(coupler_2d_bc_type), intent(in) :: IOB_struc !< Ice Ocean Boundary flux structure
-    real,                     intent(in) :: weight !< Unknown
+    real,                     intent(in) :: weight    !< A weight for accumulating these fluxes
     type(time_type), optional,intent(in) :: model_time !< Time
   end subroutine generic_tracer_coupler_accumulate
 
   !> Calls the corresponding generic_X_update_from_source routine for each package X
   subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dtts,&
        grid_dat,model_time,nbands,max_wavelength_band,sw_pen_band,opacity_band,internal_heat,&
        frunoff,grid_ht, current_wave_stress, sosga)
-    real, dimension(ilb:,jlb:,:),   intent(in) :: Temp !< Potential temperature [deg C]
-    real, dimension(ilb:,jlb:,:),   intent(in) :: Salt !< Salinity [psu]
-    real, dimension(ilb:,jlb:,:),   intent(in) :: rho_dzt !< Unknown
-    real, dimension(ilb:,jlb:,:),   intent(in) :: dzt !< Ocean layer thickness [m]
-    real, dimension(ilb:,jlb:),     intent(in) :: hblt_depth !< Boundary layer depth
-    integer,                        intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain
-    integer,                        intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain
-    integer,                        intent(in) :: tau !< Time step index of %field
-    real,                           intent(in) :: dtts !< Unknown
-    real, dimension(ilb:,jlb:),     intent(in) :: grid_dat !< Unknown
+    real, dimension(ilb:,jlb:,:),   intent(in) :: Temp   !< Potential temperature [deg C]
+    real, dimension(ilb:,jlb:,:),   intent(in) :: Salt   !< Salinity [psu]
+    real, dimension(ilb:,jlb:,:),   intent(in) :: rho_dzt !< Mass per unit area of each layer [kg m-2]
+    real, dimension(ilb:,jlb:,:),   intent(in) :: dzt    !< Ocean layer thickness [m]
+    real, dimension(ilb:,jlb:),     intent(in) :: hblt_depth !< Boundary layer depth [m]
+    integer,                        intent(in) :: ilb    !< Lower bounds of x extent of input arrays on data domain
+    integer,                        intent(in) :: jlb    !< Lower bounds of y extent of input arrays on data domain
+    integer,                        intent(in) :: tau    !< Time step index of %field
+    real,                           intent(in) :: dtts   !< The time step for this call [s]
+    real, dimension(ilb:,jlb:),     intent(in) :: grid_dat !< Grid cell areas [m2]
     type(time_type),                intent(in) :: model_time !< Time
-    integer,                        intent(in) :: nbands !< Unknown
-    real, dimension(:),             intent(in) :: max_wavelength_band !< Unknown
-    real, dimension(:,ilb:,jlb:),   intent(in) :: sw_pen_band !< Shortwave penetration
-    real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band !< Unknown
-    real, dimension(ilb:,jlb:),optional,  intent(in) :: internal_heat !< Unknown
-    real, dimension(ilb:,jlb:),optional,  intent(in) :: frunoff !< Unknown
-    real, dimension(ilb:,jlb:),optional,  intent(in) :: grid_ht !< Unknown
-    real, dimension(ilb:,jlb:),optional , intent(in) :: current_wave_stress !< Unknown
-    real,                      optional , intent(in) :: sosga !< Global average sea surface salinity
+    integer,                        intent(in) :: nbands !< The number of bands of penetrating shortwave radiation
+    real, dimension(:),             intent(in) :: max_wavelength_band !< The maximum wavelength in each band
+                                                         !! of penetrating shortwave radiation [nm]
+    real, dimension(:,ilb:,jlb:),   intent(in) :: sw_pen_band !< Penetrating shortwave radiation per band [W m-2].
+                                                         !! The wavelength or angular direction band is the first index.
+    real, dimension(:,ilb:,jlb:,:), intent(in) :: opacity_band !< Opacity of seawater averaged over each band [m-1].
+                                                         !! The wavelength or angular direction band is the first index.
+    real, dimension(ilb:,jlb:),optional,  intent(in) :: internal_heat !< Any internal or geothermal heat
+                                                         !! sources that are applied to the ocean integrated
+                                                         !! over this timestep [degC kg m-2]
+    real, dimension(ilb:,jlb:),optional,  intent(in) :: frunoff !< Rate of iceberg calving [kg m-2 s-1]
+    real, dimension(ilb:,jlb:),optional,  intent(in) :: grid_ht !< Unknown, and presently unused by MOM6
+    real, dimension(ilb:,jlb:),optional , intent(in) :: current_wave_stress !< Unknown, and presently unused by MOM6
+    real,                      optional , intent(in) :: sosga !< Global average sea surface salinity [ppt]
   end subroutine generic_tracer_source
 
   !> Update the tracers from bottom fluxes
   subroutine generic_tracer_update_from_bottom(dt, tau, model_time)
-    real,            intent(in) :: dt !< Time step increment
+    real,            intent(in) :: dt !< Time step increment [s]
     integer,         intent(in) :: tau !< Time step index used for the concentration field
     type(time_type), intent(in) :: model_time !< Time
   end subroutine generic_tracer_update_from_bottom
 
   !> Vertically diffuse all generic tracers for GOLD ocean
   subroutine generic_tracer_vertdiff_G(h_old, ea, eb, dt, kg_m2_to_H, m_to_H, tau)
-    real, dimension(:,:,:), intent(in) :: h_old !< Unknown
-    real, dimension(:,:,:), intent(in) :: ea !< Unknown
-    real, dimension(:,:,:), intent(in) :: eb !< Unknown
-    real,                   intent(in) :: dt !< Unknown
-    real,                   intent(in) :: kg_m2_to_H !< Unknown
-    real,                   intent(in) :: m_to_H !< Unknown
-    integer,                intent(in) :: tau !< Unknown
+    real, dimension(:,:,:), intent(in) :: h_old  !< Layer thickness before entrainment [H ~> m or kg m-2]
+    real, dimension(:,:,:), intent(in) :: ea     !< The amount of fluid entrained from the layer
+                                                 !! above during this call [H ~> m or kg m-2]
+    real, dimension(:,:,:), intent(in) :: eb     !< The amount of fluid entrained from the layer
+                                                 !! below during this call [H ~> m or kg m-2]
+    real,                   intent(in) :: dt     !< The amount of time covered by this call [s]
+    real,                   intent(in) :: kg_m2_to_H !< A unit conversion factor from mass per unit
+                                                 !! area to thickness units [H m2 kg-1 ~> m3 kg-1 or 1]
+    real,                   intent(in) :: m_to_H !< A unit conversion factor from heights to
+                                                 !! thickness units [H m-1 ~> 1 or kg m-3]
+    integer,                intent(in) :: tau    !< The time level to work on (always 1 for MOM6)
   end subroutine generic_tracer_vertdiff_G
 
   !> Set the coupler values for each generic tracer
@@ -115,11 +124,11 @@ subroutine generic_tracer_coupler_set(IOB_struc, ST,SS,rho,ilb,jlb,tau, dzt, sos
     integer,                     intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain
     integer,                     intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain
     integer,                     intent(in) :: tau !< Time step index of %field
-    real, dimension(ilb:,jlb:),  intent(in) :: ST !< Sea surface temperature [deg C]
-    real, dimension(ilb:,jlb:),  intent(in) :: SS !< Sea surface salinity [psu]
+    real, dimension(ilb:,jlb:),  intent(in) :: ST !< Sea surface temperature [degC]
+    real, dimension(ilb:,jlb:),  intent(in) :: SS !< Sea surface salinity [ppt]
     real, dimension(ilb:,jlb:,:,:), intent(in) :: rho !< Ocean density [kg m-3]
     real, dimension(ilb:,jlb:,:), optional, intent(in) :: dzt !< Layer thickness [m]
-    real,           optional, intent(in) :: sosga !< Unknown
+    real,           optional, intent(in) :: sosga !< Global mean sea surface salinity [ppt]
     type(time_type),optional, intent(in) :: model_time !< Time
   end subroutine generic_tracer_coupler_set