Merge pull request #229 from gustavo-marques/merge_main_08Dec22

(*) Merge GFDL to main (2022-10-27)

.github/actions/macos-setup/action.yml

@@ -13,5 +13,6 @@ runs:
         brew update
         brew install automake
         brew install netcdf
+        brew install netcdf-fortran
         brew install mpich
         echo "::endgroup::"

.github/actions/ubuntu-setup/action.yml

@@ -13,7 +13,7 @@ runs:
         sudo apt-get install netcdf-bin
         sudo apt-get install libnetcdf-dev
         sudo apt-get install libnetcdff-dev
-        sudo apt-get install mpich
-        sudo apt-get install libmpich-dev
+        sudo apt-get install openmpi-bin
+        sudo apt-get install libopenmpi-dev
         sudo apt-get install linux-tools-common
         echo "::endgroup::"

.testing/Makefile

@@ -182,14 +182,25 @@ endif
 SOURCE = \
   $(foreach ext,F90 inc c h,$(wildcard $(1)/*/*.$(ext) $(1)/*/*/*.$(ext)))
 
-MOM_SOURCE = $(call SOURCE,../src) \
-  $(wildcard ../config_src/infra/FMS1/*.F90) \
+MOM_SOURCE = \
+  $(call SOURCE,../src) \
   $(wildcard ../config_src/drivers/solo_driver/*.F90) \
   $(wildcard ../config_src/ext*/*/*.F90)
-TARGET_SOURCE = $(call SOURCE,build/target_codebase/src) \
-  $(wildcard build/target_codebase/config_src/infra/FMS1/*.F90) \
+
+TARGET_SOURCE = \
+  $(call SOURCE,build/target_codebase/src) \
   $(wildcard build/target_codebase/config_src/drivers/solo_driver/*.F90) \
   $(wildcard build/target_codebase/config_src/ext*/*.F90)
+
+# NOTE: Current default framework is FMS1, but this could change.
+ifeq ($(FRAMEWORK), fms2)
+  MOM_SOURCE +=$(wildcard ../config_src/infra/FMS2/*.F90)
+  TARGET_SOURCE += $(wildcard build/target_codebase/config_src/infra/FMS2/*.F90)
+else
+  MOM_SOURCE += $(wildcard ../config_src/infra/FMS1/*.F90)
+  TARGET_SOURCE += $(wildcard build/target_codebase/config_src/infra/FMS1/*.F90)
+endif
+
 FMS_SOURCE = $(call SOURCE,deps/fms/src)
 
 
@@ -602,6 +613,7 @@ report.cov: run.cov codecov
 	  2> build/cov/codecov.err \
 	  && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \
 	  || { \
+	    cat build/cov/codecov.err ; \
 	    echo -e "${RED}Failed to upload report.${RESET}" ; \
 		if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
 	  }
@@ -740,6 +752,7 @@ report.cov.unit: build/unit/MOM_file_parser_tests.F90.gcov codecov
 	    2> build/unit/codecov.err \
 	  && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \
 	  || { \
+	    cat build/unit/codecov.err ; \
 	    echo -e "${RED}Failed to upload report.${RESET}" ; \
 		if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
 	  }

ac/configure.ac

@@ -234,21 +234,28 @@ AC_SUBST([SRC_DIRS],
 AC_CONFIG_COMMANDS(Makefile.dep, [make depend])
 
 
-# setjmp verification
+# POSIX verification tests
 AC_LANG_PUSH([C])
 
-# Verify that either sigsetjmp (POSIX) or __sigsetjmp (glibc) are available.
-AC_CHECK_FUNC([sigsetjmp])
-AS_IF([test "$ac_cv_func_sigsetjmp" == "yes"], [
-  SIGSETJMP_NAME="sigsetjmp"
-], [
-  AC_CHECK_FUNC([__sigsetjmp], [
-    SIGSETJMP_NAME="__sigsetjmp"
-  ], [
-    AC_MSG_ERROR([Could not find a symbol for sigsetjmp.])
+# These symbols may be defined as macros, making them inaccessible by Fortran.
+# The following exist in BSD and Linux, so we just test for them.
+AC_CHECK_FUNC([setjmp], [], [AC_MSG_ERROR([Could not find setjmp.])])
+AC_CHECK_FUNC([longjmp], [], [AC_MSG_ERROR([Could not find longjmp.])])
+AC_CHECK_FUNC([siglongjmp], [], [AC_MSG_ERROR([Could not find siglongjmp.])])
+
+# Determine the sigsetjmp symbol.  If missing, then point to sigsetjmp_missing.
+#
+# Supported symbols:
+#   sigsetjmp       POSIX, BSD libc (MacOS)
+#   __sigsetjmp     glibc (Linux)
+SIGSETJMP="sigsetjmp_missing"
+for sigsetjmp_fn in sigsetjmp __sigsetjmp; do
+  AC_CHECK_FUNC([${sigsetjmp_fn}], [
+    SIGSETJMP=${sigsetjmp_fn}
+    break
   ])
-])
-AC_DEFINE_UNQUOTED([SIGSETJMP_NAME], ["$SIGSETJMP_NAME"])
+done
+AC_DEFINE_UNQUOTED([SIGSETJMP_NAME], ["${SIGSETJMP}"])
 
 # Determine the size of jmp_buf and sigjmp_buf
 AC_CHECK_SIZEOF([jmp_buf], [], [#include <setjmp.h>])

config_src/drivers/FMS_cap/ocean_model_MOM.F90

@@ -53,6 +53,7 @@ module ocean_model_mod
 use MOM_variables, only : surface
 use MOM_verticalGrid, only : verticalGrid_type
 use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
+use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
 use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
 use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
 use MOM_wave_interface, only: Update_Surface_Waves
@@ -274,9 +275,13 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
   if (.not.OS%is_ocean_pe) return
 
   OS%Time = Time_in ; OS%Time_dyn = Time_in
+  ! Call initialize MOM with an optional Ice Shelf CS which, if present triggers
+  ! initialization of ice shelf parameters and arrays.
+
   call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, &
                       OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, &
-                      diag_ptr=OS%diag, count_calls=.true., waves_CSp=OS%Waves)
+                      diag_ptr=OS%diag, count_calls=.true., ice_shelf_CSp=OS%ice_shelf_CSp, &
+                      waves_CSp=OS%Waves)
   call get_MOM_state_elements(OS%MOM_CSp, G=OS%grid, GV=OS%GV, US=OS%US, C_p=OS%C_p, &
                               C_p_scaled=OS%fluxes%C_p, use_temp=use_temperature)
 
@@ -372,9 +377,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
   endif
 
   if (OS%use_ice_shelf)  then
-    call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, &
-                              OS%diag, OS%forces, OS%fluxes)
+    call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes)
+    call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces)
   endif
+
   if (OS%icebergs_alter_ocean)  then
     call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp)
     if (.not. OS%use_ice_shelf) &

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -98,8 +98,6 @@ module MOM_cap_mod
 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
 
 public SetServices
@@ -149,7 +147,6 @@ module MOM_cap_mod
 integer              :: scalar_field_count = 0
 integer              :: scalar_field_idx_grid_nx = 0
 integer              :: scalar_field_idx_grid_ny = 0
-integer              :: nthrds  !< number of openmp threads per task
 character(len=*),parameter :: u_FILE_u = &
      __FILE__
 
@@ -465,30 +462,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   CALL ESMF_TimeIntervalGet(TINT, S=DT_OCEAN, RC=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-  !---------------------------------
-  ! openmp threads
-  !---------------------------------
-
-  call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc)
-  if (ChkErr(rc,__LINE__,u_FILE_u)) return
-
-  if (localPeCount == 1) then
-    call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, &
-         isPresent=isPresent, isSet=isSet, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    if (isPresent .and. isSet) then
-      read(cvalue,*) nthrds
-    else
-      nthrds = localPeCount
-    endif
-  else
-    nthrds = localPeCount
-  endif
-  write(logmsg,*) nthrds
-  call ESMF_LogWrite(trim(subname)//': nthreads = '//trim(logmsg), ESMF_LOGMSG_INFO)
-
-!$  call omp_set_num_threads(nthrds)
-
   call fms_init(mpi_comm_mom)
   call constants_init
   call field_manager_init
@@ -936,28 +909,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
   call ESMF_VMGet(vm, petCount=npet, mpiCommunicator=mpicom, localPet=localPet, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-  !---------------------------------
-  ! openmp threads
-  !---------------------------------
-
-  call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc)
-  if (ChkErr(rc,__LINE__,u_FILE_u)) return
-
-  if (localPeCount == 1) then
-    call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, &
-         isPresent=isPresent, isSet=isSet, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    if (isPresent .and. isSet) then
-      read(cvalue,*) nthrds
-    else
-      nthrds = localPeCount
-    endif
-  else
-    nthrds = localPeCount
-  endif
-
-!$  call omp_set_num_threads(nthrds)
-
   !---------------------------------
   ! global mom grid size
   !---------------------------------
@@ -1570,8 +1521,6 @@ subroutine ModelAdvance(gcomp, rc)
 
   call shr_file_setLogUnit (logunit)
 
-!$  call omp_set_num_threads(nthrds)
-
   ! query the Component for its clock, importState and exportState
   call ESMF_GridCompGet(gcomp, clock=clock, importState=importState, &
     exportState=exportState, rc=rc)

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -106,7 +106,7 @@ module MOM_surface_forcing
   real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [R L Z T-1 ~> Pa], if WIND_CONFIG=='gyres'
   real :: gyres_taux_n_pis   !< The number of sine lobes in the basin if WIND_CONFIG=='gyres'
   integer :: answer_date     !< This 8-digit integer gives the approximate date with which the order
-                             !! of arithmetic and and expressions were added to the code.
+                             !! of arithmetic and expressions were added to the code.
                              !! Dates before 20190101 use original answers.
                              !! Dates after 20190101 use a form of the gyre wind stresses that are
                              !! rotationally invariant and more likely to be the same between compilers.
@@ -161,8 +161,8 @@ module MOM_surface_forcing
   character(len=200) :: salinityrestore_file = '' !< The file from which to read the sea surface
                                                   !! salinity to restore toward
 
-  character(len=80)  :: stress_x_var  = '' !< X-windstress variable name in the input file
-  character(len=80)  :: stress_y_var  = '' !< Y-windstress variable name in the input file
+  character(len=80)  :: stress_x_var  = '' !< X-wind stress variable name in the input file
+  character(len=80)  :: stress_y_var  = '' !< Y-wind stress variable name in the input file
   character(len=80)  :: ustar_var     = '' !< ustar variable name in the input file
   character(len=80)  :: LW_var        = '' !< longwave heat flux variable name in the input file
   character(len=80)  :: SW_var        = '' !< shortwave heat flux variable name in the input file
@@ -447,6 +447,8 @@ subroutine wind_forcing_2gyre(sfc_state, forces, day, G, US, CS)
     forces%tauy(i,J) = 0.0
   enddo ; enddo
 
+  if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS)
+
   call callTree_leave("wind_forcing_2gyre")
 end subroutine wind_forcing_2gyre
 
@@ -484,6 +486,8 @@ subroutine wind_forcing_1gyre(sfc_state, forces, day, G, US, CS)
     forces%tauy(i,J) = 0.0
   enddo ; enddo
 
+  if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS)
+
   call callTree_leave("wind_forcing_1gyre")
 end subroutine wind_forcing_1gyre
 
@@ -499,8 +503,6 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
                                                   !! a previous surface_forcing_init call
   ! Local variables
   real :: PI            ! A common irrational number, 3.1415926535... [nondim]
-  real :: I_rho         ! The inverse of the reference density times a ratio of scaling
-                        ! factors [Z L-1 R-1 ~> m3 kg-1]
   real :: y             ! The latitude relative to the south normalized by the domain extent [nondim]
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
 
@@ -530,12 +532,7 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS)
                         forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/CS%Rho0) )
     enddo ; enddo
   else
-    I_rho = US%L_to_Z / CS%Rho0
-    do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( (CS%gust_const + &
-            sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + &
-                      (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho )
-    enddo ; enddo
+    call stresses_to_ustar(forces, G, US, CS)
   endif
 
   call callTree_leave("wind_forcing_gyres")
@@ -558,8 +555,6 @@ subroutine Neverworld_wind_forcing(sfc_state, forces, day, G, US, CS)
   real :: Pa_to_RLZ_T2  ! A unit conversion factor from Pa to the internal units
                         ! for wind stresses [R Z L T-2 Pa-1 ~> 1]
   real :: PI            ! A common irrational number, 3.1415926535... [nondim]
-  real :: I_rho         ! The inverse of the reference density times a ratio of scaling
-                        ! factors [Z L-1 R-1 ~> m3 kg-1]
   real :: y             ! The latitude relative to the south normalized by the domain extent [nondim]
   real :: tau_max       ! The magnitude of the wind stress [R Z L T-2 ~> Pa]
   real :: off           ! An offset in the relative latitude [nondim]
@@ -602,14 +597,7 @@ subroutine Neverworld_wind_forcing(sfc_state, forces, day, G, US, CS)
   enddo ; enddo
 
   ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1].  ustar is always positive.
-  if (associated(forces%ustar)) then
-    I_rho = US%L_to_Z / CS%Rho0
-    do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = sqrt( (CS%gust_const + &
-            sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + &
-                      (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho )
-    enddo ; enddo
-  endif
+  if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS)
 
 end subroutine Neverworld_wind_forcing
 
@@ -625,8 +613,6 @@ subroutine scurve_wind_forcing(sfc_state, forces, day, G, US, CS)
                                                     !! a previous surface_forcing_init call
   ! Local variables
   integer :: i, j, kseg
-  real :: I_rho         ! The inverse of the reference density times a ratio of scaling
-                        ! factors [Z L-1 R-1 ~> m3 kg-1]
   real :: y_curve       ! The latitude relative to the southern end of a curve segment [degreesN]
   real :: L_curve       ! The latitudinal extent of a curve segment [degreesN]
 ! real :: ydata(7) = (/ -70., -45., -15., 0., 15., 45., 70. /)
@@ -657,14 +643,7 @@ subroutine scurve_wind_forcing(sfc_state, forces, day, G, US, CS)
   enddo ; enddo
 
   ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1].  ustar is always positive.
-  if (associated(forces%ustar)) then
-    I_rho = US%L_to_Z / CS%Rho0
-    do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      forces%ustar(i,j) = sqrt( (CS%gust_const + &
-            sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + &
-                      (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho )
-    enddo ; enddo
-  endif
+  if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS)
 
 end subroutine scurve_wind_forcing
 
@@ -892,6 +871,37 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS)
   call callTree_leave("wind_forcing_by_data_override")
 end subroutine wind_forcing_by_data_override
 
+!> Translate the wind stresses into the friction velocity, including effects of background gustiness.
+subroutine stresses_to_ustar(forces, G, US, CS)
+  type(mech_forcing),       intent(inout) :: forces !< A structure with the driving mechanical forces
+  type(ocean_grid_type),    intent(in)    :: G      !< Grid structure.
+  type(unit_scale_type),    intent(in)    :: US     !< A dimensional unit scaling type
+  type(surface_forcing_CS), pointer       :: CS     !< pointer to control structure returned by
+                                                    !! a previous surface_forcing_init call
+  ! Local variables
+  real :: I_rho         ! The inverse of the reference density times a ratio of scaling
+                        ! factors [Z L-1 R-1 ~> m3 kg-1]
+  integer :: i, j, is, ie, js, je
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+
+  I_rho = US%L_to_Z / CS%Rho0
+
+  if (CS%read_gust_2d) then
+    do j=js,je ; do i=is,ie
+      forces%ustar(i,j) = sqrt( (CS%gust(i,j) + &
+              sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + &
+                        (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * I_rho )
+    enddo ; enddo
+  else
+    do j=js,je ; do i=is,ie
+      forces%ustar(i,j) = sqrt( (CS%gust_const + &
+            sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + &
+                      (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho )
+    enddo ; enddo
+  endif
+
+end subroutine stresses_to_ustar
 
 !> Specifies zero surface buoyancy fluxes from input files.
 subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS)