Merge branch 'dev/ncar' into restructure_diabatic_driver

.gitlab-ci.yml

@@ -116,8 +116,9 @@ run:
     - time tar zxf $CACHE_DIR/build-pgi-repro-$CI_PIPELINE_ID.tgz
     # time tar zxf $CACHE_DIR/build-gnu-debug-$CI_PIPELINE_ID.tgz
     - echo "make -f MRS/Makefile.tests all -B" > job.sh
-    - msub -l partition=c4,nodes=29,walltime=00:29:00,qos=norm -q debug -S /bin/tcsh -j oe -A gfdl_o -z -o log.$CI_PIPELINE_ID -N mom6_regression -K job.sh
+    - msub -l partition=c4,nodes=29,walltime=00:31:00,qos=norm -q debug -S /bin/tcsh -j oe -A gfdl_o -z -o log.$CI_PIPELINE_ID -N mom6_regression -K job.sh
     - cat log.$CI_PIPELINE_ID
+    - test -f restart_results_gnu.tar.gz
     - time tar zvcf $CACHE_DIR/results-$CI_PIPELINE_ID.tgz *.tar.gz
 
 # Tests

config_src/coupled_driver/MOM_surface_forcing.F90

@@ -492,6 +492,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, &
         fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
       enddo ; enddo
     endif
+    fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
   endif
 
   ! more salt restoring logic
@@ -645,6 +646,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
     call safe_alloc_ptr(forces%rigidity_ice_v,isd,ied,JsdB,JedB)
   endif
 
+  forces%accumulate_rigidity = .true. ! Multiple components may contribute to rigidity.
   if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
   if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
 
@@ -661,7 +663,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
         forces%p_surf(i,j) = forces%p_surf_full(i,j)
       enddo ; enddo
     endif
+  else
+    do j=js,je ; do i=is,ie
+      forces%p_surf_full(i,j) = 0.0
+      forces%p_surf(i,j) = 0.0
+    enddo ; enddo
   endif
+  forces%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
 
   wind_stagger = CS%wind_stagger
   if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. &

config_src/coupled_driver/ocean_model_MOM.F90

@@ -54,7 +54,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 : ice_shelf_end, ice_shelf_save_restart
+use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
 use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type
 use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums
 use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data
@@ -514,18 +514,24 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
                              OS%grid, OS%forcing_CSp)
 
   if (OS%fluxes%fluxes_used) then
-    call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, &
+    if (do_thermo) &
+      call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, &
                                OS%grid, OS%forcing_CSp, OS%sfc_state, &
                                OS%restore_salinity, OS%restore_temp)
 
     ! Add ice shelf fluxes
     if (OS%use_ice_shelf) then
-      call shelf_calc_flux(OS%sfc_state, OS%forces, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
+      if (do_thermo) &
+        call shelf_calc_flux(OS%sfc_state, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
+      if (do_dyn) &
+        call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces)
     endif
     if (OS%icebergs_alter_ocean)  then
-      call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, &
-                          OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
-      call iceberg_fluxes(OS%grid, OS%fluxes, OS%use_ice_shelf, &
+      if (do_dyn) &
+        call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, &
+                            OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
+      if (do_thermo) &
+        call iceberg_fluxes(OS%grid, OS%fluxes, OS%use_ice_shelf, &
                           OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
     endif
 
@@ -541,22 +547,28 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     OS%fluxes%dt_buoy_accum = dt_coupling
   else
     OS%flux_tmp%C_p = OS%fluxes%C_p
-    call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, &
+    if (do_thermo) &
+      call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, &
                                OS%grid, OS%forcing_CSp, OS%sfc_state, OS%restore_salinity,OS%restore_temp)
 
     if (OS%use_ice_shelf) then
-      call shelf_calc_flux(OS%sfc_state, OS%forces, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
+      if (do_thermo) &
+        call shelf_calc_flux(OS%sfc_state, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
+      if (do_dyn) &
+        call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces)
     endif
     if (OS%icebergs_alter_ocean)  then
-      call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, &
-                          OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
-      call iceberg_fluxes(OS%grid, OS%flux_tmp, OS%use_ice_shelf, &
+      if (do_dyn) &
+        call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, &
+                            OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
+      if (do_thermo) &
+        call iceberg_fluxes(OS%grid, OS%flux_tmp, OS%use_ice_shelf, &
                           OS%sfc_state, dt_coupling, OS%marine_ice_CSp)
     endif
 
     call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, dt_coupling, OS%grid, weight)
     ! Some of the fields that exist in both the forcing and mech_forcing types
-    ! are time-averages must be copied back to the forces type.
+    ! (e.g., ustar) are time-averages must be copied back to the forces type.
     call copy_back_forcing_fields(OS%fluxes, OS%forces, OS%grid)
 
 #ifdef _USE_GENERIC_TRACER

config_src/mct_driver/ocn_comp_mct.F90

@@ -61,7 +61,7 @@ module ocn_comp_mct
 use MOM_diag_mediator,    only: diag_mediator_close_registration, diag_mediator_end
 use MOM_diag_mediator,    only: safe_alloc_ptr
 use MOM_ice_shelf,        only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
-use MOM_ice_shelf,        only: ice_shelf_end, ice_shelf_save_restart
+use MOM_ice_shelf,        only: add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
 use MOM_string_functions, only: uppercase
 use MOM_constants,        only: CELSIUS_KELVIN_OFFSET, hlf, hlv
 use MOM_EOS,              only: gsw_sp_from_sr, gsw_pt_from_ct
@@ -1727,7 +1727,8 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, &
 
     ! Add ice shelf fluxes
     if (OS%use_ice_shelf) then
-      call shelf_calc_flux(OS%sfc_state, OS%forces, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp)
+      call shelf_calc_flux(OS%sfc_state, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp)
+      call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces)
     endif
 
     ! GMM, check ocean_model_MOM.F90 to enable the following option
@@ -1748,7 +1749,8 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, &
                        OS%restore_salinity,OS%restore_temp)
 
     if (OS%use_ice_shelf) then
-      call shelf_calc_flux(OS%sfc_state, OS%forces, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp)
+      call shelf_calc_flux(OS%sfc_state, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp)
+      call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces)
     endif
 
     ! GMM, check ocean_model_MOM.F90 to enable the following option
@@ -1947,6 +1949,7 @@ subroutine ocn_import(forces, fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp,
     fluxes%heat_added(:,:)=0.0
     fluxes%salt_flux_added(:,:)=0.0
   endif
+  forces%accumulate_rigidity = .true. ! Multiple components may contribute to rigidity.
   if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0
   if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0
 
@@ -2142,6 +2145,7 @@ subroutine ocn_import(forces, fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp,
       else
         forces%p_surf_SSH => forces%p_surf_full
       endif
+      forces%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
 
     endif
 

config_src/solo_driver/MOM_driver.F90

@@ -66,7 +66,7 @@ program MOM_main
   use time_interp_external_mod, only : time_interp_external_init
 
   use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS
-  use MOM_ice_shelf, only : shelf_calc_flux, ice_shelf_save_restart
+  use MOM_ice_shelf, only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart
 ! , add_shelf_flux_forcing, add_shelf_flux_IOB
 
   use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
@@ -246,8 +246,8 @@ program MOM_main
   endif
 
 !$  call omp_set_num_threads(ocean_nthreads)
-!$OMP PARALLEL private(adder)
 !$  base_cpu = get_cpu_affinity()
+!$OMP PARALLEL private(adder)
 !$  if (use_hyper_thread) then
 !$     if (mod(omp_get_thread_num(),2) == 0) then
 !$        adder = omp_get_thread_num()/2
@@ -258,7 +258,7 @@ program MOM_main
 !$     adder = omp_get_thread_num()
 !$  endif
 !$  call set_cpu_affinity (base_cpu + adder)
-!$  write(6,*) " ocean  ", omp_get_num_threads(), get_cpu_affinity(), adder, omp_get_thread_num()
+!$  write(6,*) " ocean ", base_cpu, get_cpu_affinity(), adder, omp_get_thread_num(), omp_get_num_threads()
 !$  call flush(6)
 !$OMP END PARALLEL
 
@@ -483,10 +483,8 @@ program MOM_main
     endif
 
     if (use_ice_shelf) then
-      call shelf_calc_flux(sfc_state, forces, fluxes, Time, dt_forcing, ice_shelf_CSp)
-!###IS     call add_shelf_flux_forcing(fluxes, ice_shelf_CSp)
-!###IS  ! With a coupled ice/ocean run, use the following call.
-!###IS      call add_shelf_flux_IOB(ice_ocean_bdry_type, ice_shelf_CSp)
+      call shelf_calc_flux(sfc_state, fluxes, Time, dt_forcing, ice_shelf_CSp)
+      call add_shelf_forces(grid, Ice_shelf_CSp, forces)
     endif
     fluxes%fluxes_used = .false.
     fluxes%dt_buoy_accum = dt_forcing

src/core/MOM.F90

@@ -52,9 +52,9 @@ module MOM
 use MOM_ALE,                   only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile
 use MOM_ALE,                   only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
-use MOM_legacy_diabatic_driver,only : legacy_diabatic
 use MOM_diabatic_driver,       only : diabatic, diabatic_driver_init, diabatic_CS
 use MOM_diabatic_driver,       only : adiabatic, adiabatic_driver_init, diabatic_driver_end
+use MOM_diabatic_driver,       only : legacy_diabatic
 use MOM_diagnostics,           only : calculate_diagnostic_fields, MOM_diagnostics_init
 use MOM_diagnostics,           only : register_transport_diags, post_transport_diagnostics
 use MOM_diagnostics,           only : register_surface_diags, write_static_fields
@@ -2994,8 +2994,7 @@ subroutine MOM_end(CS)
   call tracer_registry_end(CS%tracer_Reg)
   call tracer_flow_control_end(CS%tracer_flow_CSp)
 
-  ! GMM, the following is commented because it fails on Travis.
-  !if (associated(CS%diabatic_CSp)) call diabatic_driver_end(CS%diabatic_CSp)
+  call diabatic_driver_end(CS%diabatic_CSp)
 
   if (CS%offline_tracer_mode) call offline_transport_end(CS%offline_CSp)
 

src/core/MOM_forcing_type.F90

@@ -118,6 +118,10 @@ module MOM_forcing_type
                                !! in corrections to the sea surface height field
                                !! that is passed back to the calling routines.
                                !! This may point to p_surf or to p_surf_full.
+  logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere
+                                 !! and various types of ice needs to be accumulated, and the
+                                 !! surface pressure explicitly reset to zero at the driver level
+                                 !! when appropriate.
 
   ! tide related inputs
   real, pointer, dimension(:,:) :: &
@@ -207,6 +211,13 @@ module MOM_forcing_type
                                  !< enabled, and is exactly 0 away from shelves or on land.
     rigidity_ice_u => NULL(), &  !< Depth-integrated lateral viscosity of ice
     rigidity_ice_v => NULL()     !< shelves or sea ice at u- or v-points (m3/s)
+  logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere
+                                 !! and various types of ice needs to be accumulated, and the
+                                 !! surface pressure explicitly reset to zero at the driver level
+                                 !! when appropriate.
+  logical :: accumulate_rigidity = .false. !< If true, the rigidity due to various types of
+                                 !! ice needs to be accumulated, and the rigidity explicitly
+                                 !! reset to zero at the driver level when appropriate.
 
   logical :: initialized = .false. !< This indicates whether the appropriate
                                  !! arrays have been initialized.

src/core/MOM_grid.F90

@@ -574,8 +574,8 @@ end subroutine MOM_grid_end
 !! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points.
 !! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point.
 !!
-!! \image html Grid_metrics.png
-!! "The labelling of distances (grid metrics) at various staggered location on an T-cell and around a q-point.
+!! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered
+!! location on an T-cell and around a q-point."
 !!
 !! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively.
 !!

src/core/MOM_open_boundary.F90

@@ -1814,10 +1814,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt)
          J=segment%HI%JsdB
          allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
          do k=1,nz
-           rx_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k)
-           rx_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k)
+           rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k)
+           rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k)
            do I=segment%HI%IsdB+1,segment%HI%IedB-1
-             rx_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k))
+             rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k))
            enddo
          enddo
          if (segment%radiation_tan) then
@@ -1925,10 +1925,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt)
          J=segment%HI%JsdB
          allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
          do k=1,nz
-           rx_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k)
-           rx_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k)
+           rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k)
+           rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k)
            do I=segment%HI%IsdB+1,segment%HI%IedB-1
-             rx_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k))
+             rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k))
            enddo
          enddo
          if (segment%radiation_tan) then