Merge pull request mom-ocean#21 from gustavo-marques/merge_dev_gfdl_1…

…4feb2022

Merge dev/gfdl as of 14 Feb 2022

.testing/Makefile

@@ -146,6 +146,7 @@ endif
 # These are set to true by our Travis configuration if testing a pull request
 DO_REGRESSION_TESTS ?=
 REPORT_COVERAGE ?=
+CODECOV_UPLOADER_URL ?= https://uploader.codecov.io/latest/linux/codecov
 
 ifeq ($(DO_REGRESSION_TESTS), true)
   BUILDS += target
@@ -165,6 +166,7 @@ else
 endif
 
 
+
 # List of source files to link this Makefile's dependencies to model Makefiles
 # Assumes a depth of two, and the following extensions: F90 inc c h
 # (1): Root directory
@@ -542,6 +544,7 @@ $(foreach c,$(CONFIGS),$(eval $(call CONFIG_DIM_RULE,$(c))))
 # $(4): MOM_override configuration
 # $(5): Environment variables
 # $(6): Number of MPI ranks
+
 define STAT_RULE
 work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 $(VENV_PATH)
 	@echo "Running test $$*.$(1)..."
@@ -570,10 +573,13 @@ work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 $(VENV_PATH)
 	@echo -e "$(DONE): $$*.$(1); no runtime errors."
 	if [ $(3) ]; then \
 	  mkdir -p results/$$* ; \
-	  cd build/symmetric \
-	    && bash <(curl -s https://codecov.io/bash) -Z -n $$@ \
-	      > codecov.$$*.$(1).out \
-	      2> codecov.$$*.$(1).err \
+	  cd build/symmetric ; \
+	  gcov *.gcda > gcov.$$*.$(1).out ; \
+	  curl -s $(CODECOV_UPLOADER_URL) -o codecov ; \
+	  chmod +x codecov ; \
+	  ./codecov -Z -f "*.gcov" -n $$@ \
+	    > codecov.$$*.$(1).out \
+	    2> codecov.$$*.$(1).err \
 	    && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}"; \
 	fi
 endef

README.md

@@ -1,6 +1,5 @@
-[![Build Status](https://travis-ci.org/NOAA-GFDL/MOM6.svg?branch=dev/master)](https://travis-ci.org/NOAA-GFDL/MOM6)
 [![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/ocean_model_MOM.F90

@@ -86,7 +86,7 @@ module ocean_model_mod
 
 !> This type is used for communication with other components via the FMS coupler.
 !! The element names and types can be changed only with great deliberation, hence
-!! the persistnce of things like the cutsy element name "avg_kount".
+!! the persistence of things like the cutesy element name "avg_kount".
 type, public ::  ocean_public_type
   type(domain2d) :: Domain    !< The domain for the surface fields.
   logical :: is_ocean_pe      !< .true. on processors that run the ocean model.
@@ -110,8 +110,8 @@ module ocean_model_mod
                     !! a global max across ocean and non-ocean processors can be
                     !! used to determine its value.
   real, pointer, dimension(:,:)  :: &
-    t_surf => NULL(), & !< SST on t-cell (degrees Kelvin)
-    s_surf => NULL(), & !< SSS on t-cell (psu)
+    t_surf => NULL(), & !< SST on t-cell [degrees Kelvin]
+    s_surf => NULL(), & !< SSS on t-cell [ppt]
     u_surf => NULL(), & !< i-velocity at the locations indicated by stagger [m s-1].
     v_surf => NULL(), & !< j-velocity at the locations indicated by stagger [m s-1].
     sea_lev => NULL(), & !< Sea level in m after correction for surface pressure,
@@ -154,8 +154,8 @@ module ocean_model_mod
 
   logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
                               !! ocean dynamics and forcing fluxes.
-  real :: press_to_z          !< A conversion factor between pressure and ocean
-                              !! depth in m, usually 1/(rho_0*g) [m Pa-1].
+  real :: press_to_z          !< A conversion factor between pressure and ocean depth,
+                              !! usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1].
   real :: C_p                 !< The heat capacity of seawater [J degC-1 kg-1].
   logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode
                               !! with the barotropic and baroclinic dynamics, thermodynamics,
@@ -221,7 +221,7 @@ module ocean_model_mod
 !! for restarts and reading restart files if appropriate.
 !!
 !!   This subroutine initializes both the ocean state and the ocean surface type.
-!! Because of the way that indicies and domains are handled, Ocean_sfc must have
+!! Because of the way that indices and domains are handled, Ocean_sfc must have
 !! been used in a previous call to initialize_ocean_type.
 subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas_fields_ocn)
   type(ocean_public_type), target, &
@@ -242,16 +242,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
   ! Local variables
-  real :: Rho0        ! The Boussinesq ocean density [kg m-3].
-  real :: G_Earth     ! The gravitational acceleration [m s-2].
-  real :: HFrz        !< If HFrz > 0 (m), melt potential will be computed.
+  real :: Rho0        ! The Boussinesq ocean density [R ~> kg m-3]
+  real :: G_Earth     ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
+  real :: HFrz        !< If HFrz > 0 [Z ~> m], melt potential will be computed.
                       !! The actual depth over which melt potential is computed will
                       !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
                       !! If HFrz <= 0 (default), melt potential will not be computed.
-  logical :: use_melt_pot!< If true, allocate melt_potential array
+  logical :: use_melt_pot !< If true, allocate melt_potential array
 
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   character(len=40)  :: mdl = "ocean_model_init"  ! This module's name.
   character(len=48)  :: stagger ! A string indicating the staggering locations for the
                                 ! surface velocities returned to the coupler.
@@ -331,28 +331,29 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                  "calculate accelerations and the mass for conservation "//&
                  "properties, or with BOUSSINSEQ false to convert some "//&
                  "parameters from vertical units of m to kg m-2.", &
-                 units="kg m-3", default=1035.0)
+                 units="kg m-3", default=1035.0, scale=OS%US%kg_m3_to_R)
   call get_param(param_file, mdl, "G_EARTH", G_Earth, &
                  "The gravitational acceleration of the Earth.", &
-                 units="m s-2", default = 9.80)
+                 units="m s-2", default=9.80, scale=OS%US%m_s_to_L_T**2*OS%US%Z_to_m)
 
   call get_param(param_file, mdl, "ICE_SHELF",  OS%use_ice_shelf, &
                  "If true, enables the ice shelf model.", default=.false.)
 
   call get_param(param_file, mdl, "ICEBERGS_APPLY_RIGID_BOUNDARY",  OS%icebergs_alter_ocean, &
                  "If true, allows icebergs to change boundary condition felt by ocean", default=.false.)
 
-  OS%press_to_z = 1.0/(Rho0*G_Earth)
+  OS%press_to_z = 1.0 / (Rho0*G_Earth)
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
   call get_param(param_file, mdl, "HFREEZE", HFrz, &
                  "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
                  "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
                  "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
-                 "melt potential will not be computed.", units="m", default=-1.0, do_not_log=.true.)
+                 "melt potential will not be computed.", &
+                 units="m", default=-1.0, scale=OS%US%m_to_Z, do_not_log=.true.)
 
-  if (HFrz .gt. 0.0) then
+  if (HFrz > 0.0) then
     use_melt_pot=.true.
   else
     use_melt_pot=.false.
@@ -655,7 +656,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
 
 ! Translate state into Ocean.
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
-!                                   Ice_ocean_boundary%p, OS%press_to_z)
+!                                   OS%fluxes%p_surf_full, OS%press_to_z)
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
   Time1 = OS%Time ; if (do_dyn) Time1 = OS%Time_dyn
   call coupler_type_send_data(Ocean_sfc%fields, Time1)
@@ -766,7 +767,7 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_field
                                               !! tracer fluxes.
 
   integer :: xsz, ysz, layout(2)
-  ! ice-ocean-boundary fields are always allocated using absolute indicies
+  ! ice-ocean-boundary fields are always allocated using absolute indices
   ! and have no halos.
   integer :: isc, iec, jsc, jec
 
@@ -806,7 +807,7 @@ end subroutine initialize_ocean_public_type
 !! surface state variable.  This may eventually be folded into the MOM
 !! code that calculates the surface state in the first place.
 !! Note the offset in the arrays because the ocean_data_type has no
-!! halo points in its arrays and always uses absolute indicies.
+!! halo points in its arrays and always uses absolute indices.
 subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_to_z)
   type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                !! describe the surface state of the ocean.
@@ -816,9 +817,9 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
                                                !! have their data set here.
   type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
   type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
-  real,        optional, intent(in)    :: patm(:,:)  !< The pressure at the ocean surface [Pa].
-  real,        optional, intent(in)    :: press_to_z !< A conversion factor between pressure and
-                                               !! ocean depth in m, usually 1/(rho_0*g) [m Pa-1].
+  real,        optional, intent(in)    :: patm(:,:)  !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
+  real,        optional, intent(in)    :: press_to_z !< A conversion factor between pressure and ocean
+                                               !! depth, usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1]
   ! Local variables
   real :: IgR0
   character(len=48)  :: val_str
@@ -860,7 +861,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
 
   if (present(patm)) then
     do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
-      Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z
+      Ocean_sfc%sea_lev(i,j) = US%Z_to_m * (sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z)
       Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0)
     enddo ; enddo
   else
@@ -946,7 +947,7 @@ end subroutine ocean_model_init_sfc
 
 !> ocean_model_flux_init is used to initialize properties of the air-sea fluxes
 !! as determined by various run-time parameters.  It can be called from
-!! non-ocean PEs, or PEs that have not yet been initialzed, and it can safely
+!! non-ocean PEs, or PEs that have not yet been initialized, and it can safely
 !! be called multiple times.
 subroutine ocean_model_flux_init(OS, verbosity)
   type(ocean_state_type), optional, pointer :: OS  !< An optional pointer to the ocean state,

config_src/drivers/mct_cap/ocn_comp_mct.F90

@@ -722,7 +722,7 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn)
   call mct_gGrid_importRattr(dom_ocn,"lat",data,lsize)
 
   k = 0
-  L2_to_rad2 = grid%US%L_to_m**2 / grid%Rad_Earth**2
+  L2_to_rad2 = 1.0 / grid%Rad_Earth_L**2
   do j = grid%jsc, grid%jec
     do i = grid%isc, grid%iec
       k = k + 1 ! Increment position within gindex

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
@@ -1111,7 +1112,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
            k = k + 1 ! Increment position within gindex
            if (mask(k) /= 0) then
               mesh_areas(k) = dataPtr_mesh_areas(k)
-              model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
+              model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth_L**2
               mod2med_areacor(k) = model_areas(k) / mesh_areas(k)
               med2mod_areacor(k) = mesh_areas(k) / model_areas(k)
            end if
@@ -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"