Merge pull request #65 from jiandewang/feature/update-to-GFDL-20210402

Feature/update to gfdl 20210402

.testing/Makefile

@@ -84,6 +84,9 @@ FCFLAGS_COVERAGE ?=
 # - FMS cannot be built with the same aggressive initialization flags as MOM6,
 #   so FCFLAGS_INIT is used to provide additional MOM6 configuration.
 
+# User-defined LDFLAGS (applied to all builds and FMS)
+LDFLAGS_USER ?=
+
 # Set to `true` to require identical results from DEBUG and REPRO builds
 # NOTE: Many compilers (Intel, GCC on ARM64) do not yet produce identical
 #   results across DEBUG and REPRO builds (as defined below), so we disable on
@@ -217,8 +220,8 @@ REPRO_FCFLAGS := FCFLAGS="$(FCFLAGS_REPRO) $(FCFLAGS_FMS)"
 OPENMP_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)"
 TARGET_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)"
 
-MOM_LDFLAGS := LDFLAGS="$(LDFLAGS_FMS)"
-SYMMETRIC_LDFLAGS := LDFLAGS="$(COVERAGE) $(LDFLAGS_FMS)"
+MOM_LDFLAGS := LDFLAGS="$(LDFLAGS_FMS) $(LDFLAGS_USER)"
+SYMMETRIC_LDFLAGS := LDFLAGS="$(COVERAGE) $(LDFLAGS_FMS) $(LDFLAGS_USER)"
 
 
 # Environment variable configuration
@@ -286,7 +289,7 @@ $(TARGET_CODEBASE)/ac/configure: $(TARGET_CODEBASE)
 
 $(TARGET_CODEBASE):
 	git clone --recursive $(MOM_TARGET_URL) $@
-	cd $@ && git checkout $(MOM_TARGET_BRANCH)
+	cd $@ && git checkout --recurse-submodules $(MOM_TARGET_BRANCH)
 
 
 #---

config_src/drivers/ice_solo_driver/ice_shelf_driver.F90

@@ -185,6 +185,8 @@ program Shelf_main
     endif
   endif
 
+  call Get_MOM_Input(param_file, dirs)
+
   ! Read ocean_solo restart, which can override settings from the namelist.
   if (file_exists(trim(dirs%restart_input_dir)//'ice_solo.res')) then
     call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ice_solo.res', action=READONLY_FILE)
@@ -215,7 +217,6 @@ program Shelf_main
     Start_time = real_to_time(0.0)
   endif
 
-  call Get_MOM_Input(param_file, dirs)
   ! Determining the internal unit scaling factors for this run.
   call unit_scaling_init(param_file, US)
 

config_src/drivers/mct_cap/mom_surface_forcing_mct.F90

@@ -496,7 +496,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     if (associated(fluxes%frunoff)) then
       fluxes%latent(i,j)              = fluxes%latent(i,j) - &
           IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
-      fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
+      fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) &
+          * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
     endif
     ! contribution from evaporation
     if (associated(IOB%q_flux)) then

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -136,7 +136,7 @@ module MOM_cap_mod
 logical              :: profile_memory = .true.
 logical              :: grid_attach_area = .false.
 logical              :: use_coldstart = .true.
-logical              :: use_mommesh = .false.
+logical              :: use_mommesh = .true.
 character(len=128)   :: scalar_field_name = ''
 integer              :: scalar_field_count = 0
 integer              :: scalar_field_idx_grid_nx = 0
@@ -349,7 +349,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc)
   write(logmsg,*) use_coldstart
   call ESMF_LogWrite('MOM_cap:use_coldstart = '//trim(logmsg), ESMF_LOGMSG_INFO)
 
-  use_mommesh = .false.
+  use_mommesh = .true.
   call NUOPC_CompAttributeGet(gcomp, name="use_mommesh", value=value, &
        isPresent=isPresent, isSet=isSet, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return

config_src/drivers/solo_driver/MOM_driver.F90

@@ -250,6 +250,10 @@ program MOM_main
   ! This call sets the number and affinity of threads with openMP.
   !$  call set_MOM_thread_affinity(ocean_nthreads, use_hyper_thread)
 
+  ! This call is required to initiate dirs%restart_input_dir for ocean_solo.res
+  ! The contents of dirs will be reread in initialize_MOM.
+  call get_MOM_input(dirs=dirs)
+
   ! Read ocean_solo restart, which can override settings from the namelist.
   if (file_exists(trim(dirs%restart_input_dir)//'ocean_solo.res')) then
     call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ocean_solo.res', action=READONLY_FILE)

config_src/infra/FMS1/MOM_diag_manager_infra.F90

@@ -20,9 +20,9 @@ module MOM_diag_manager_infra
 use diag_manager_mod, only : register_diag_field_fms => register_diag_field
 use diag_manager_mod, only : register_static_field_fms => register_static_field
 use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id
-use time_manager_mod, only : time_type
+use MOM_time_manager, only : time_type
 use MOM_domain_infra, only : MOM_domain_type
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_error_infra,  only : MOM_error => MOM_err, FATAL, WARNING
 
 implicit none ; private
 

config_src/infra/FMS1/MOM_domain_infra.F90

@@ -4,7 +4,7 @@ module MOM_domain_infra
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_coms_infra,  only : PE_here, root_PE, num_PEs
-use MOM_cpu_clock,   only : cpu_clock_begin, cpu_clock_end
+use MOM_cpu_clock_infra, only : cpu_clock_begin, cpu_clock_end
 use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, WARNING, FATAL
 
 use mpp_domains_mod, only : domain2D, domain1D
@@ -1689,6 +1689,8 @@ subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, &
   if ((MD_in%io_layout(1) + MD_in%io_layout(2) > 0) .and. &
       (MD_in%layout(1)*MD_in%layout(2) > 1)) then
     call mpp_define_io_domain(mpp_domain, MD_in%io_layout)
+  else
+    call mpp_define_io_domain(mpp_domain, (/ 1, 1 /) )
   endif
 
 end subroutine clone_MD_to_d2D

src/core/MOM_isopycnal_slopes.F90

@@ -37,14 +37,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
                                                                      !! thermodynamic variables
   real,                                        intent(in)    :: dt_kappa_smooth !< A smoothing vertical diffusivity
                                                                      !! times a smoothing timescale [Z2 ~> m2].
-  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim]
-  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim]
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
                                      optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at
-                                                                     !! interfaces between u-points [T-2 ~> s-2]
+                                                                     !! interfaces between u-points [L2 Z-2 T-2 ~> s-2]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
                                      optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at
-                                                                     !! interfaces between u-points [T-2 ~> s-2]
+                                                                     !! interfaces between v-points [L2 Z-2 T-2 ~> s-2]
   integer,                           optional, intent(in)    :: halo !< Halo width over which to compute
   type(ocean_OBC_type),              optional, pointer       :: OBC  !< Open boundaries control structure.
 
@@ -86,15 +86,15 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
   real :: drdz          ! Vertical density gradient [R Z-1 ~> kg m-4].
   real :: Slope         ! The slope of density surfaces, calculated in a way
                         ! that is always between -1 and 1.
-  real :: mag_grad2     ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
+  real :: mag_grad2     ! The squared magnitude of the 3-d density gradient [R2 Z-2 ~> kg2 m-8].
   real :: slope2_Ratio  ! The ratio of the slope squared to slope_max squared.
   real :: h_neglect     ! A thickness that is so small it is usually lost
                         ! in roundoff and can be neglected [H ~> m or kg m-2].
   real :: h_neglect2    ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
   real :: dz_neglect    ! A change in interface heighs that is so small it is usually lost
                         ! in roundoff and can be neglected [Z ~> m].
   logical :: use_EOS    ! If true, density is calculated from T & S using an equation of state.
-  real :: G_Rho0        ! The gravitational acceleration divided by density [Z2 T-2 R-1 ~> m5 kg-2 s-2]
+  real :: G_Rho0        ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
   real :: Z_to_L        ! A conversion factor between from units for e to the
                         ! units for lateral distances.
   real :: L_to_Z        ! A conversion factor between from units for lateral distances
@@ -134,7 +134,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
 
   present_N2_u = PRESENT(N2_u)
   present_N2_v = PRESENT(N2_v)
-  G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / GV%Rho0
+  G_Rho0 = GV%g_Earth / GV%Rho0
   if (present_N2_u) then
     do j=js,je ; do I=is-1,ie
       N2_u(I,j,1) = 0.
@@ -248,17 +248,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
 
         ! This estimate of slope is accurate for small slopes, but bounded
         ! to be between -1 and 1.
-        mag_grad2 = drdx**2 + (L_to_Z*drdz)**2
+        mag_grad2 = (Z_to_L*drdx)**2 + drdz**2
         if (mag_grad2 > 0.0) then
           slope_x(I,j,K) = drdx / sqrt(mag_grad2)
         else ! Just in case mag_grad2 = 0 ever.
           slope_x(I,j,K) = 0.0
         endif
 
-        if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy frequency [T-2 ~> s-2]
+        if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
 
       else ! With .not.use_EOS, the layers are constant density.
-        slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j)
+        slope_x(I,j,K) = (e(i,j,K)-e(i+1,j,K)) * G%IdxCu(I,j)
       endif
       if (local_open_u_BC) then
         l_seg = OBC%segnum_u(I,j)
@@ -351,17 +351,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
 
         ! This estimate of slope is accurate for small slopes, but bounded
         ! to be between -1 and 1.
-        mag_grad2 = drdy**2 + (L_to_Z*drdz)**2
+        mag_grad2 = (Z_to_L*drdy)**2 + drdz**2
         if (mag_grad2 > 0.0) then
           slope_y(i,J,K) = drdy / sqrt(mag_grad2)
         else ! Just in case mag_grad2 = 0 ever.
           slope_y(i,J,K) = 0.0
         endif
 
-        if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy frequency [T-2 ~> s-2]
+        if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
 
       else ! With .not.use_EOS, the layers are constant density.
-        slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J)
+        slope_y(i,J,K) = (e(i,j,K)-e(i,j+1,K)) * G%IdyCv(i,J)
       endif
       if (local_open_v_BC) then
         l_seg = OBC%segnum_v(i,J)

src/ice_shelf/MOM_ice_shelf.F90

@@ -722,6 +722,17 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
     endif
   endif
 
+  ! Melting has been computed, now is time to update thickness and mass with dynamic ice shelf
+  if (CS%active_shelf_dynamics) then
+    call change_thickness_using_melt(ISS, G, US, US%s_to_T*time_step, fluxes, CS%density_ice, CS%debug)
+
+    if (CS%debug) then
+      call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m)
+      call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, &
+                   scale=US%RZ_to_kg_m2)
+    endif
+  endif
+
   if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0)
 
   call add_shelf_flux(G, US, CS, sfc_state, fluxes)