merged to head of dev/ncar

.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

@@ -45,9 +45,7 @@ module MOM_surface_forcing
 
 #include <MOM_memory.h>
 
-public IOB_allocate
-public convert_IOB_to_fluxes
-public convert_IOB_to_forces
+public convert_IOB_to_fluxes, convert_IOB_to_forces
 public surface_forcing_init
 public ice_ocn_bnd_type_chksum
 public forcing_save_restart
@@ -494,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
@@ -647,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
 
@@ -663,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. &
@@ -841,58 +847,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
   call cpu_clock_end(id_clock_forcing)
 end subroutine convert_IOB_to_forces
 
-subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
-
-  type(ice_ocean_boundary_type), intent(inout)    :: IOB    !< An ice-ocean boundary type with fluxes to drive
-  integer, intent(in) :: isc, iec, jsc, jec                 !< The ocean's local grid size
-
-    allocate ( IOB% u_flux (isc:iec,jsc:jec),          &
-               IOB% v_flux (isc:iec,jsc:jec),          &
-               IOB% t_flux (isc:iec,jsc:jec),          &
-               IOB% q_flux (isc:iec,jsc:jec),          &
-               IOB% salt_flux (isc:iec,jsc:jec),       &
-               IOB% lw_flux (isc:iec,jsc:jec),         &
-               IOB% sw_flux_vis_dir (isc:iec,jsc:jec), &
-               IOB% sw_flux_vis_dif (isc:iec,jsc:jec), &
-               IOB% sw_flux_nir_dir (isc:iec,jsc:jec), &
-               IOB% sw_flux_nir_dif (isc:iec,jsc:jec), &
-               IOB% lprec (isc:iec,jsc:jec),           &
-               IOB% fprec (isc:iec,jsc:jec),           &
-               IOB% runoff (isc:iec,jsc:jec),          &
-               IOB% ustar_berg (isc:iec,jsc:jec),      &
-               IOB% area_berg (isc:iec,jsc:jec),       &
-               IOB% mass_berg (isc:iec,jsc:jec),       &
-               IOB% calving (isc:iec,jsc:jec),         &
-               IOB% runoff_hflx (isc:iec,jsc:jec),     &
-               IOB% calving_hflx (isc:iec,jsc:jec),    &
-               IOB% mi (isc:iec,jsc:jec),              &
-               IOB% p (isc:iec,jsc:jec))
-
-    IOB%u_flux          = 0.0
-    IOB%v_flux          = 0.0
-    IOB%t_flux          = 0.0
-    IOB%q_flux          = 0.0
-    IOB%salt_flux       = 0.0
-    IOB%lw_flux         = 0.0
-    IOB%sw_flux_vis_dir = 0.0
-    IOB%sw_flux_vis_dif = 0.0
-    IOB%sw_flux_nir_dir = 0.0
-    IOB%sw_flux_nir_dif = 0.0
-    IOB%lprec           = 0.0
-    IOB%fprec           = 0.0
-    IOB%runoff          = 0.0
-    IOB%ustar_berg      = 0.0
-    IOB%area_berg       = 0.0
-    IOB%mass_berg       = 0.0
-    IOB%calving         = 0.0
-    IOB%runoff_hflx     = 0.0
-    IOB%calving_hflx    = 0.0
-    IOB%mi              = 0.0
-    IOB%p               = 0.0
-
-end subroutine IOB_allocate
-
-
 !> Adds thermodynamic flux adjustments obtained via data_override
 !! Component name is 'OCN'
 !! Available adjustments are:

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
@@ -77,7 +77,6 @@ module ocean_model_mod
 implicit none ; private
 
 public ocean_model_init, ocean_model_end, update_ocean_model
-public get_ocean_grid ! add by Jiande
 public ocean_model_save_restart, Ocean_stock_pe
 public ice_ocean_boundary_type
 public ocean_model_init_sfc, ocean_model_flux_init
@@ -515,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
 
@@ -542,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
@@ -1093,28 +1104,8 @@ subroutine ocean_model_data2D_get(OS,Ocean, name, array2D,isc,jsc)
      array2D(isc:,jsc:) = Ocean%t_surf(isc:,jsc:)-CELSIUS_KELVIN_OFFSET
   case('btfHeat')
      array2D(isc:,jsc:) = 0
-  case('tlat')
-     array2D(isc:,jsc:) = OS%grid%geoLatT(g_isc:g_iec,g_jsc:g_jec)
-  case('tlon')
-     array2D(isc:,jsc:) = OS%grid%geoLonT(g_isc:g_iec,g_jsc:g_jec)
-  case('ulat')
-     array2D(isc:,jsc:) = OS%grid%geoLatCu(g_isc:g_iec,g_jsc:g_jec)
-  case('ulon')
-     array2D(isc:,jsc:) = OS%grid%geoLonCu(g_isc:g_iec,g_jsc:g_jec)
-  case('vlat')
-     array2D(isc:,jsc:) = OS%grid%geoLatCv(g_isc:g_iec,g_jsc:g_jec)
-  case('vlon')
-     array2D(isc:,jsc:) = OS%grid%geoLonCv(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLatBu')
-     array2D(isc:,jsc:) = OS%grid%geoLatBu(g_isc:g_iec,g_jsc:g_jec)
-  case('geoLonBu')
-     array2D(isc:,jsc:) = OS%grid%geoLonBu(g_isc:g_iec,g_jsc:g_jec)
-  case('cos_rot')
-     array2D(isc:,jsc:) = OS%grid%cos_rot(g_isc:g_iec,g_jsc:g_jec) ! =1
-  case('sin_rot')
-     array2D(isc:,jsc:) = OS%grid%sin_rot(g_isc:g_iec,g_jsc:g_jec) ! =0
   case default
-     call MOM_error(FATAL,'ocean_model_data2D_get: unknown argument name='//name)
+     call MOM_error(FATAL,'get_ocean_grid_data2D: unknown argument name='//name)
   end select
 
 
@@ -1133,7 +1124,7 @@ subroutine ocean_model_data1D_get(OS,Ocean, name, value)
   case('c_p')
      value = OS%C_p
   case default
-     call MOM_error(FATAL,'ocean_model_data1D_get: unknown argument name='//name)
+     call MOM_error(FATAL,'get_ocean_grid_data1D: unknown argument name='//name)
   end select
 
 
@@ -1161,21 +1152,4 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
 
 end subroutine ocean_public_type_chksum
 
-!#######################################################################
-! <SUBROUTINE NAME="get_ocean_grid">
-!
-! <DESCRIPTION>
-! Obtain the ocean grid.
-! </DESCRIPTION>
-!
-  subroutine get_ocean_grid(OS, Gridp)
-    type(ocean_state_type) :: OS
-    type(ocean_grid_type) , pointer          :: Gridp
-
-    Gridp => OS%grid
-    return
-
-  end subroutine get_ocean_grid
-! </SUBROUTINE> NAME="get_ocean_grid"
-
 end module ocean_model_mod

config_src/mct_driver/ocn_comp_mct.F90

@@ -8,7 +8,8 @@ module ocn_comp_mct
 use ESMF,                only: ESMF_ClockGet, ESMF_TimeGet, ESMF_TimeIntervalGet
 use seq_cdata_mod,       only: seq_cdata, seq_cdata_setptrs
 use seq_flds_mod,        only: seq_flds_x2o_fields, seq_flds_o2x_fields 
-use shr_flds_mod,        only: shr_flds_dom_coord, shr_flds_dom_other
+!use seq_flds_mod,       only: shr_flds_dom_coord=>seq_flds_dom_coord, shr_flds_dom_other=>seq_flds_dom_other
+use seq_flds_mod,       only: shr_flds_dom_coord, shr_flds_dom_other
 use mct_mod,             only: mct_gsMap, mct_gsmap_init, mct_gsMap_lsize, &
                                mct_gsmap_orderedpoints
 use mct_mod,             only: mct_aVect, mct_aVect_init, mct_aVect_zero, &
@@ -490,7 +491,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
 
   ! Translate import fields to ice_ocean_boundary
   !TODO: make this an input variable
-  glb%sw_decomp = .false.
+  !glb%sw_decomp = .false.
   !END TODO:
   if (glb%sw_decomp) then
      call ocn_import(x2o_o%rattr, glb%ind,  glb%grid, Ice_ocean_boundary, glb%ocn_public, glb%stdout, Eclock, &
@@ -520,7 +521,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
     write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') trim(runid), year, month, day, seconds
 
     call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, &
-                      glb%ocn_state%restart_CSp, .false., filename=restartname,GV=glb%ocn_state%GV)
+                      glb%ocn_state%restart_CSp, .false., filename=restartname, GV=glb%ocn_state%GV)
 
     ! write name of restart file in the rpointer file
     nu = shr_file_getUnit()

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

@@ -54,6 +54,7 @@ module MOM
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
 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
@@ -197,6 +198,9 @@ module MOM
 
   logical :: adiabatic               !< If true, there are no diapycnal mass fluxes, and no calls
                                      !! to routines to calculate or apply diapycnal fluxes.
+  logical :: use_legacy_diabatic_driver!< If true (default), use the a legacy version of the
+                                       !! diabatic subroutine. This is temporary and is needed
+                                       !! to avoid change in answers.
   logical :: diabatic_first          !< If true, apply diabatic and thermodynamic
                                      !! processes before time stepping the dynamics.
   logical :: use_ALE_algorithm       !< If true, use the ALE algorithm rather than layered
@@ -1151,8 +1155,14 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, &
     endif
 
     call cpu_clock_begin(id_clock_diabatic)
-    call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
-                  dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    if (CS%use_legacy_diabatic_driver) then
+      ! the following subroutine is legacy and will be deleted in the near future.
+      call legacy_diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
+                          dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    else
+      call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
+                    dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    endif
     fluxes%fluxes_used = .true.
     call cpu_clock_end(id_clock_diabatic)
 
@@ -1627,6 +1637,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "true. This assumes that KD = KDML = 0.0 and that \n"//&
                  "there is no buoyancy forcing, but makes the model \n"//&
                  "faster by eliminating subroutine calls.", default=.false.)
+  call get_param(param_file, "MOM", "USE_LEGACY_DIABATIC_DRIVER", CS%use_legacy_diabatic_driver, &
+                 "If true, use the a legacy version of the diabatic subroutine. \n"//&
+                 "This is temporary and is needed avoid change in answers.", &
+                 default=.true.)
   call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, &
                  "If False, skips the dynamics calls that update u & v, as well as \n"//&
                  "the gravity wave adjustment to h. This is a fragile feature and \n"//&
@@ -2364,6 +2378,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   if (associated(CS%visc%Kv_shear)) &
     call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1)
 
+  if (associated(CS%visc%Kv_slow)) &
+    call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1)
+
   call cpu_clock_end(id_clock_pass_init)
 
   call register_obsolete_diagnostics(param_file, CS%diag)
@@ -2977,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)