Merge branch 'dev/gfdl' into internal_tides_split_freqs_modes

.github/workflows/coverage.yml

@@ -28,10 +28,14 @@ jobs:
     - name: Report unit test coverage to CI (PR)
       if: github.event_name == 'pull_request'
       run: make report.cov.unit REQUIRE_COVERAGE_UPLOAD=true
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
 
     - name: Report unit test coverage to CI (Push)
       if: github.event_name != 'pull_request'
       run: make report.cov.unit
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
 
     - name: Compile ocean-only MOM6 with code coverage
       run: make -j build/cov/MOM6
@@ -42,7 +46,11 @@ jobs:
     - name: Report coverage to CI (PR)
       if: github.event_name == 'pull_request'
       run: make report.cov REQUIRE_COVERAGE_UPLOAD=true
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
 
     - name: Report coverage to CI (Push)
       if: github.event_name != 'pull_request'
       run: make report.cov
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

.testing/Makefile

@@ -571,13 +571,21 @@ endef
 
 # Upload coverage reports
 CODECOV_UPLOADER_URL ?= https://uploader.codecov.io/latest/linux/codecov
+CODECOV_TOKEN ?=
+
+ifdef CODECOV_TOKEN
+  CODECOV_TOKEN_ARG = -t $(CODECOV_TOKEN)
+else
+  CODECOV_TOKEN_ARG =
+endif
+
 codecov:
 	curl -s $(CODECOV_UPLOADER_URL) -o $@
 	chmod +x codecov
 
 .PHONY: report.cov
 report.cov: run.cov codecov
-	./codecov -R build/cov -Z -f "*.gcov" \
+	./codecov $(CODECOV_TOKEN_ARG) -R build/cov -Z -f "*.gcov" \
 	  > build/cov/codecov.out \
 	  2> build/cov/codecov.err \
 	  && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \
@@ -690,7 +698,7 @@ build/unit/MOM_file_parser_tests.F90.gcov: $(WORKSPACE)/work/unit/std.out
 
 .PHONY: report.cov.unit
 report.cov.unit: build/unit/MOM_file_parser_tests.F90.gcov codecov
-	./codecov -R build/unit -f "*.gcov" -Z -n "Unit tests" \
+	./codecov $(CODECOV_TOKEN_ARG) -R build/unit -f "*.gcov" -Z -n "Unit tests" \
 	    > build/unit/codecov.out \
 	    2> build/unit/codecov.err \
 	  && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \

ac/makedep

@@ -150,7 +150,10 @@ def create_deps(src_dirs, makefile, debug, exec_target, fc_rule,
                 dep for pair in zip(found_mods, found_objs) for dep in pair
             ]
             missing_mods = [m for m in o2uses[o] if m not in all_modules]
-            incs = nested_inc(o2h[o] + o2inc[o], f2F)
+
+            incs, inc_used = nested_inc(o2h[o] + o2inc[o], f2F)
+            inc_mods = [u for u in inc_used if u not in found_mods and u in all_modules]
+
             incdeps = sorted(set([f2F[f] for f in incs if f in f2F]))
             incargs = sorted(set(['-I'+os.path.dirname(f) for f in incdeps]))
             if debug:
@@ -167,7 +170,7 @@ def create_deps(src_dirs, makefile, debug, exec_target, fc_rule,
                 print("#   program:", ' '.join(o2prg[o]), file=file)
             if o2mods[o]:
                 print(' '.join(o2mods[o])+':', o, file=file)
-            print(o + ':', o2F90[o], ' '.join(incdeps+found_deps), file=file)
+            print(o + ':', o2F90[o], ' '.join(inc_mods + incdeps + found_deps), file=file)
             print('\t'+fc_rule, ' '.join(incargs), file=file)
 
         # Write rule for each object from C
@@ -243,21 +246,30 @@ def link_obj(obj, o2uses, mod2o, all_modules):
 def nested_inc(inc_files, f2F):
     """List of all files included by "inc_files", either by #include or F90
     include."""
+    hlst = []
+    used_mods = set()
+
     def recur(hfile):
         if hfile not in f2F.keys():
             return
-        _, _, cpp, inc, _, _ = scan_fortran_file(f2F[hfile])
+
+        _, used, cpp, inc, _, _ = scan_fortran_file(f2F[hfile])
+
+        # Record any module updates inside of include files
+        used_mods.update(used)
+
         if len(cpp) + len(inc) > 0:
             for h in cpp+inc:
                 if h not in hlst and h in f2F.keys():
                     recur(h)
                     hlst.append(h)
             return
         return
-    hlst = []
+
     for h in inc_files:
         recur(h)
-    return inc_files + sorted(set(hlst))
+
+    return inc_files + sorted(set(hlst)), used_mods
 
 
 def scan_fortran_file(src_file):
@@ -268,8 +280,10 @@ def scan_fortran_file(src_file):
         lines = file.readlines()
 
         external_namespace = True
+            # True if we are in the external (i.e. global) namespace
 
         file_has_externals = False
+            # True if the file contains any external objects
 
         for line in lines:
             match = re_module.match(line.lower())
@@ -321,17 +335,18 @@ def object_file(src_file):
 def find_files(src_dirs):
     """Return sorted list of all source files starting from each directory in
     the list "src_dirs"."""
+
+    # TODO: Make this a user-defined argument
+    extensions = ('.f90', '.f', '.c', '.inc', '.h',  '.fh')
+
     files = []
+
     for path in src_dirs:
         if not os.path.isdir(path):
             raise ValueError("Directory '{}' was not found".format(path))
         for p, d, f in os.walk(os.path.normpath(path), followlinks=True):
             for file in f:
-                # TODO: use any()
-                if (file.endswith('.F90') or file.endswith('.f90')
-                        or file.endswith('.f') or file.endswith('.F')
-                        or file.endswith('.h') or file.endswith('.inc')
-                        or file.endswith('.c') or file.endswith('.H')):
+                if any(file.lower().endswith(ext) for ext in extensions):
                     files.append(p+'/'+file)
     return sorted(set(files))
 

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -63,6 +63,7 @@ module MOM_surface_forcing_gfdl
                                 !! the winds that are being provided in calls to
                                 !! update_ocean_model.
   logical :: use_temperature    !< If true, temp and saln used as state variables.
+  logical :: nonBous            !< If true, this run is fully non-Boussinesq
   real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim].
 
   real :: Rho0                  !< Boussinesq reference density [R ~> kg m-3]
@@ -282,8 +283,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   ! allocation and initialization if this is the first time that this
   ! flux type has been used.
   if (fluxes%dt_buoy_accum < 0) then
-    call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., press=.true., &
-                               fix_accum_bug=CS%fix_ustar_gustless_bug, tau_mag=.true.)
+    call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.not.CS%nonBous, press=.true., &
+                               fix_accum_bug=CS%fix_ustar_gustless_bug, tau_mag=CS%nonBous)
 
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -718,8 +719,8 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   ! allocation and initialization if this is the first time that this
   ! mechanical forcing type has been used.
   if (.not.forces%initialized) then
-    call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., &
-                               press=.true., tau_mag=.true.)
+    call allocate_mech_forcing(G, forces, stress=.true., ustar=.not.CS%nonBous, &
+                               press=.true., tau_mag=CS%nonBous)
 
     call safe_alloc_ptr(forces%p_surf,isd,ied,jsd,jed)
     call safe_alloc_ptr(forces%p_surf_full,isd,ied,jsd,jed)
@@ -792,14 +793,26 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   ! Set the wind stresses and ustar.
   if (wt1 <= 0.0) then
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, &
-                              ustar=forces%ustar, mag_tau=forces%tau_mag, tau_halo=1)
+                              tau_halo=1)
+    if (associated(forces%ustar)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=forces%ustar)
+    if (associated(forces%tau_mag)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, mag_tau=forces%tau_mag)
   else
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux=forces%taux, tauy=forces%tauy, &
-                              ustar=ustar_tmp, mag_tau=tau_mag_tmp, tau_halo=1)
-    do j=js,je ; do i=is,ie
-      forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
-      forces%tau_mag(i,j) = wt1*forces%tau_mag(i,j) + wt2*tau_mag_tmp(i,j)
-    enddo ; enddo
+                              tau_halo=1)
+    if (associated(forces%ustar)) then
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=ustar_tmp)
+      do j=js,je ; do i=is,ie
+        forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
+      enddo ; enddo
+    endif
+    if (associated(forces%tau_mag)) then
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, mag_tau=tau_mag_tmp)
+      do j=js,je ; do i=is,ie
+        forces%tau_mag(i,j) = wt1*forces%tau_mag(i,j) + wt2*tau_mag_tmp(i,j)
+      enddo ; enddo
+    endif
   endif
 
   ! Find the net mass source in the input forcing without other adjustments.
@@ -960,7 +973,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
 
   ! Set surface momentum stress related fields as a function of staggering.
   if (present(taux) .or. present(tauy) .or. &
-      ((do_ustar.or.do_gustless) .and. .not.associated(IOB%stress_mag)) ) then
+      ((do_ustar .or. do_tau_mag .or. do_gustless) .and. .not.associated(IOB%stress_mag)) ) then
 
     if (wind_stagger == BGRID_NE) then
       taux_in_B(:,:) = 0.0 ; tauy_in_B(:,:) = 0.0
@@ -1278,6 +1291,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
   ! Local variables
   real :: utide             ! The RMS tidal velocity [Z T-1 ~> m s-1].
   real :: Flux_const_dflt   ! A default piston velocity for restoring surface properties [m day-1]
+  logical :: Boussinesq       ! If true, this run is fully Boussinesq
+  logical :: semi_Boussinesq  ! If true, this run is partially non-Boussinesq
   real :: rho_TKE_tidal     ! The constant bottom density used to translate tidal amplitudes into the
                             ! tidal bottom TKE input used with INT_TIDE_DISSIPATION [R ~> kg m-3]
   logical :: new_sim              ! False if this simulation was started from a restart file
@@ -1320,12 +1335,20 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
   call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, &
                  "If true, Temperature and salinity are used as state "//&
                  "variables.", default=.true.)
+  call get_param(param_file, mdl, "BOUSSINESQ", Boussinesq, &
+                 "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
+  call get_param(param_file, mdl, "SEMI_BOUSSINESQ", semi_Boussinesq, &
+                 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
+                 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
+                 "height changes.  This only applies if BOUSSINESQ is false.", &
+                 default=.true., do_not_log=.true.)
+  CS%nonBous = .not.(Boussinesq .or. semi_Boussinesq)
   call get_param(param_file, mdl, "RHO_0", CS%Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to "//&
                  "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, scale=US%kg_m3_to_R)
+                 units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) ! (, do_not_log=CS%nonBous)
   call get_param(param_file, mdl, "LATENT_HEAT_FUSION", CS%latent_heat_fusion, &
                  "The latent heat of fusion.", units="J/kg", default=hlf, scale=US%J_kg_to_Q)
   call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", CS%latent_heat_vapor, &

config_src/drivers/ice_solo_driver/ice_shelf_driver.F90

@@ -24,7 +24,7 @@ program Shelf_main
   use MOM_cpu_clock,       only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
   use MOM_cpu_clock,       only : CLOCK_COMPONENT
   use MOM_debugging,       only : MOM_debugging_init
-  use MOM_diag_mediator,   only : diag_mediator_init, diag_mediator_infrastructure_init
+  use MOM_diag_mediator,   only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info
   use MOM_diag_mediator,   only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration
   use MOM_domains,         only : MOM_infra_init, MOM_infra_end
   use MOM_domains,         only : MOM_domains_init, clone_MOM_domain, pass_var
@@ -54,6 +54,8 @@ program Shelf_main
   use MOM_verticalGrid,    only : verticalGrid_type, verticalGridInit, verticalGridEnd
   use MOM_write_cputime,   only : write_cputime, MOM_write_cputime_init
   use MOM_write_cputime,   only : write_cputime_start_clock, write_cputime_CS
+  use MOM_forcing_type,    only : forcing
+  use MOM_ice_shelf_initialize, only : initialize_ice_SMB
 
   use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS
   use MOM_ice_shelf, only : ice_shelf_save_restart, solo_step_ice_shelf
@@ -75,7 +77,9 @@ program Shelf_main
   ! CPU time limit.  nmax is determined by evaluating the CPU time used between successive calls to
   ! write_cputime.  Initially it is set to be very large.
   integer :: nmax=2000000000
-
+  ! A structure containing pointers to the thermodynamic forcing fields
+  ! at the ocean surface.
+  type(forcing) :: fluxes
   ! A structure containing several relevant directory paths.
   type(directories) :: dirs
 
@@ -104,7 +108,7 @@ program Shelf_main
   real :: time_step               ! The time step [T ~> s]
 
   ! A pointer to a structure containing metrics and related information.
-  type(ocean_grid_type), pointer :: ocn_grid
+  type(ocean_grid_type), pointer :: ocn_grid => NULL()
 
   type(dyn_horgrid_type), pointer :: dG => NULL()   ! A dynamic version of the horizontal grid
   type(hor_index_type),   pointer :: HI => NULL()   ! A hor_index_type for array extents
@@ -114,7 +118,7 @@ program Shelf_main
   type(ocean_OBC_type),          pointer :: OBC => NULL()
 
   ! A pointer to a structure containing dimensional unit scaling factors.
-  type(unit_scale_type), pointer :: US
+  type(unit_scale_type), pointer :: US => NULL()
 
   type(diag_ctrl), pointer :: &
     diag => NULL()              ! A pointer to the diagnostic regulatory structure
@@ -138,8 +142,9 @@ program Shelf_main
   integer :: yr, mon, day, hr, mins, sec   ! Temp variables for writing the date.
   type(param_file_type) :: param_file      ! The structure indicating the file(s)
                                            ! containing all run-time parameters.
+  real :: smb !A constant surface mass balance that can be specified in the param_file
   character(len=9)  :: month
-  character(len=16) :: calendar = 'julian'
+  character(len=16) :: calendar = 'noleap'
   integer :: calendar_type=-1
 
   integer :: unit, io_status, ierr
@@ -184,6 +189,8 @@ program Shelf_main
     endif
   endif
 
+  ! Get the names of the I/O directories and initialization file.
+  ! Also calls the subroutine that opens run-time parameter files.
   call Get_MOM_Input(param_file, dirs)
 
   ! Read ocean_solo restart, which can override settings from the namelist.
@@ -252,8 +259,11 @@ program Shelf_main
   ! Set up the ocean model domain and grid; the ice model grid is set in initialize_ice_shelf,
   ! but the grids have strong commonalities in this configuration, and the ocean grid is required
   ! to set up the diag mediator control structure.
-  call MOM_domains_init(ocn_grid%domain, param_file)
+  allocate(ocn_grid)
+  call MOM_domains_init(ocn_grid%domain, param_file) !, domain_name='MOM')
+  allocate(HI)
   call hor_index_init(ocn_grid%Domain, HI, param_file)
+  allocate(dG)
   call create_dyn_horgrid(dG, HI)
   call clone_MOM_domain(ocn_grid%Domain, dG%Domain)
 
@@ -266,11 +276,16 @@ program Shelf_main
   ! Initialize the diag mediator.  The ocean's vertical grid is not really used here, but at
   ! present the interface to diag_mediator_init assumes the presence of ocean-specific information.
   call verticalGridInit(param_file, GV, US)
+  allocate(diag)
   call diag_mediator_init(ocn_grid, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
 
   call callTree_waypoint("returned from diag_mediator_init()")
 
-  call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag)
+  call set_axes_info(ocn_grid, GV, US, param_file, diag)
+
+  call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag, fluxes_in=fluxes, solo_ice_sheet_in=.true.)
+
+  call initialize_ice_SMB(fluxes%shelf_sfc_mass_flux, ocn_grid, US, param_file)
 
   ! This is the end of the code that is the counterpart of MOM_initialization.
   call callTree_waypoint("End of ice shelf initialization.")
@@ -378,7 +393,7 @@ program Shelf_main
 
     ! This call steps the model over a time time_step.
     Time1 = Master_Time ; Time = Master_Time
-    call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time)
+    call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time, fluxes_in=fluxes)
 
 !   Time = Time + Time_step_shelf
 !   This is here to enable fractional-second time steps.
@@ -412,6 +427,20 @@ program Shelf_main
       if (BTEST(Restart_control,0)) then
         call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir)
       endif
+      ! Write ice shelf solo restart file.
+      if (is_root_pe())then
+        call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res')
+        write(unit, '(i6,8x,a)') calendar_type, &
+          '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)'
+
+        call get_date(Start_time, yr, mon, day, hr, mins, sec)
+        write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, &
+          'Model start time:   year, month, day, hour, minute, second'
+        call get_date(Time, yr, mon, day, hr, mins, sec)
+        write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, &
+          'Current model time: year, month, day, hour, minute, second'
+        call close_file(unit)
+      endif
       restart_time = restart_time + restint
     endif
 
@@ -456,12 +485,11 @@ program Shelf_main
   endif
 
   call callTree_waypoint("End Shelf_main")
+  call ice_shelf_end(ice_shelf_CSp)
   call diag_mediator_end(Time, diag, end_diag_manager=.true.)
   if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.)
   call cpu_clock_end(termClock)
 
   call io_infra_end ; call MOM_infra_end
 
-  call ice_shelf_end(ice_shelf_CSp)
-
 end program Shelf_main