Tc4 ice shelf

.testing/tc4/build_grid.py

@@ -2,14 +2,13 @@
 from netCDF4 import stringtochar
 import numpy as np
 
-nx, ny = 14, 10     # Grid size
+nx, ny = 10, 8     # Grid size
 depth0 = 100.       # Uniform depth
 ds = 0.01           # grid resolution at the equator in degrees
 Re = 6.378e6        # Radius of earth
 
 topo_ = np.zeros((ny, nx)) + depth0
 f_topo = nc.Dataset('topog.nc', 'w', format='NETCDF3_CLASSIC')
-ny, nx = topo_.shape
 f_topo.createDimension('ny', ny)
 f_topo.createDimension('nx', nx)
 f_topo.createDimension('ntiles', 1)
@@ -26,8 +25,9 @@
 dx = np.zeros((2*ny + 1, 2*nx))
 dy = np.zeros((2*ny, 2*nx + 1))
 rad_deg = np.pi / 180.
-dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1])
-         * np.cos(0.5*rad_deg*(y[:, 0:-1] + y[:, 1:])))
+#dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1])
+#         * np.cos(0.5*rad_deg*(y[:, 0:-1] + y[:, 1:])))
+dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1]))
 dy[:] = rad_deg * Re * (y[1:, :] - y[0:-1, :])
 
 f_sg = nc.Dataset('ocean_hgrid.nc', 'w', format='NETCDF3_CLASSIC')
@@ -74,3 +74,17 @@
 f_sg.variables['tile'][:] = str_
 f_sg.sync()
 f_sg.close()
+
+thick_ = np.zeros((ny, nx))
+area_  = area[::2,::2]+area[::2,1::2]+area[1::2,::2]+area[1::2,1::2]
+thick_[:,:int(nx/2)]=0.5*depth0
+area_[thick_==0.]=0.
+f_thick = nc.Dataset('shelf.nc', 'w', format='NETCDF3_CLASSIC')
+f_thick.createDimension('ny', ny)
+f_thick.createDimension('nx', nx)
+f_thick.createVariable('thick', 'f8', ('ny', 'nx'))
+f_thick.createVariable('area', 'f8', ('ny', 'nx'))
+f_thick.variables['thick'][:] = thick_
+f_thick.variables['area'][:] = area_
+f_thick.sync()
+f_thick.close()

config_src/solo_driver/MOM_driver.F90

@@ -28,6 +28,7 @@ program MOM_main
   use MOM_cpu_clock,       only : CLOCK_COMPONENT
   use MOM_diag_mediator,   only : enable_averaging, disable_averaging, diag_mediator_end
   use MOM_diag_mediator,   only : diag_ctrl, diag_mediator_close_registration
+  use MOM_IS_diag_mediator,   only : diag_IS_ctrl=>diag_ctrl, diag_mediator_IS_end=>diag_mediator_end
   use MOM,                 only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
   use MOM,                 only : extract_surface_state, finish_MOM_initialization
   use MOM,                 only : get_MOM_state_elements, MOM_state_is_synchronized
@@ -61,7 +62,7 @@ program MOM_main
   use MOM_verticalGrid,    only : verticalGrid_type
   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_get_input,       only : get_MOM_input
   use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_size
   use ensemble_manager_mod, only : ensemble_pelist_setup
   use mpp_mod, only : set_current_pelist => mpp_set_current_pelist
@@ -80,22 +81,22 @@ program MOM_main
 #include <MOM_memory.h>
 
   ! A structure with the driving mechanical surface forces
-  type(mech_forcing) :: forces
+  type(mech_forcing), pointer :: forces => NULL()
   ! A structure containing pointers to the thermodynamic forcing fields
   ! at the ocean surface.
-  type(forcing) :: fluxes
+  type(forcing), pointer :: fluxes => NULL()
 
   ! A structure containing pointers to the ocean surface state fields.
-  type(surface) :: sfc_state
+  type(surface), pointer :: sfc_state => NULL()
 
   ! A pointer to a structure containing metrics and related information.
-  type(ocean_grid_type), pointer :: grid
-  type(verticalGrid_type), pointer :: GV
+  type(ocean_grid_type), pointer :: grid => NULL()
+  type(verticalGrid_type), pointer :: GV => NULL()
   ! A pointer to a structure containing dimensional unit scaling factors.
-  type(unit_scale_type), pointer :: US
+  type(unit_scale_type), pointer :: US => NULL()
 
   ! If .true., use the ice shelf model for part of the domain.
-  logical :: use_ice_shelf
+  logical :: use_ice_shelf = .false.
 
   ! If .true., use surface wave coupling
   logical :: use_waves = .false.
@@ -198,8 +199,10 @@ program MOM_main
   type(MOM_restart_CS),      pointer :: &
     restart_CSp => NULL()     !< A pointer to the restart control structure
                               !! that will be used for MOM restart files.
-  type(diag_ctrl), pointer :: &
-    diag => NULL()            !< A pointer to the diagnostic regulatory structure
+  type(diag_ctrl),           pointer :: &
+       diag => NULL()         !< A pointer to the diagnostic regulatory structure
+  type(diag_IS_ctrl), pointer :: &
+      diag_IS => NULL()       !< A pointer to the diagnostic regulatory structure
   !-----------------------------------------------------------------------
 
   character(len=4), parameter :: vers_num = 'v2.0'
@@ -219,6 +222,8 @@ program MOM_main
 
   call MOM_infra_init() ; call io_infra_init()
 
+  allocate(forces,fluxes,sfc_state)
+
   ! Initialize the ensemble manager.  If there are no settings for ensemble_size
   ! in input.nml(ensemble.nml), these should not do anything.  In coupled
   ! configurations, this all occurs in the external driver.
@@ -297,16 +302,34 @@ program MOM_main
     ! In this case, the segment starts at a time fixed by ocean_solo.res
     segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6))
     Time = segment_start_time
-    call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, &
-                        segment_start_time, offline_tracer_mode=offline_tracer_mode, &
-                        diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp)
   else
     ! In this case, the segment starts at a time read from the MOM restart file
     ! or left as Start_time by MOM_initialize.
     Time = Start_time
+  endif
+
+  ! Read paths and filenames from namelist and store in "dirs".
+  ! Also open the parsed input parameter file(s) and setup param_file.
+  call get_MOM_input(param_file, dirs)
+
+  call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, &
+       "If true, enables the ice shelf model.", default=.false.)
+  if (use_ice_shelf) then
+    ! These arrays are not initialized in most solo cases, but are needed
+    ! when using an ice shelf
+    call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, &
+                              diag_IS, forces, fluxes, sfc_state)
+  endif
+  call close_param_file(param_file)
+
+  if (sum(date) >= 0) then
+    call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, &
+                        segment_start_time, offline_tracer_mode=offline_tracer_mode, &
+                        diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp)
+  else
     call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, &
                         offline_tracer_mode=offline_tracer_mode, diag_ptr=diag, &
-                        tracer_flow_CSp=tracer_flow_CSp)
+                        tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp)
   endif
 
   call get_MOM_state_elements(MOM_CSp, G=grid, GV=GV, US=US, C_p_scaled=fluxes%C_p)
@@ -320,14 +343,6 @@ program MOM_main
                             surface_forcing_CSp, tracer_flow_CSp)
   call callTree_waypoint("done surface_forcing_init")
 
-  call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, &
-                 "If true, enables the ice shelf model.", default=.false.)
-  if (use_ice_shelf) then
-    ! These arrays are not initialized in most solo cases, but are needed
-    ! when using an ice shelf
-    call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, &
-                              diag, forces, fluxes)
-  endif
 
   call get_param(param_file,mod_name,"USE_WAVES",Use_Waves,&
        "If true, enables surface wave modules.",default=.false.)
@@ -434,6 +449,7 @@ program MOM_main
   call close_param_file(param_file)
   call diag_mediator_close_registration(diag)
 
+
   ! Write out a time stamp file.
   if (calendar_type /= NO_CALENDAR) then
     call open_file(unit, 'time_stamp.out', form=ASCII_FILE, action=APPEND_FILE, &
@@ -482,7 +498,7 @@ program MOM_main
 
     if (use_ice_shelf) then
       call shelf_calc_flux(sfc_state, fluxes, Time, dt_forcing, ice_shelf_CSp)
-      call add_shelf_forces(grid, US, Ice_shelf_CSp, forces)
+      call add_shelf_forces(US, Ice_shelf_CSp, forces, external_call=.true.)
     endif
     fluxes%fluxes_used = .false.
     fluxes%dt_buoy_accum = US%s_to_T*dt_forcing
@@ -652,6 +668,9 @@ program MOM_main
 
   call callTree_waypoint("End MOM_main")
   call diag_mediator_end(Time, diag, end_diag_manager=.true.)
+  if (use_ice_shelf) then
+     call diag_mediator_IS_end(Time, diag_IS)
+  endif
   if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.)
   call cpu_clock_end(termClock)
 

src/ALE/MOM_ALE.F90

@@ -331,7 +331,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
   type(ALE_CS),                               pointer       :: CS  !< Regridding parameters and options
   type(ocean_OBC_type),                       pointer       :: OBC !< Open boundary structure
   real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main [T ~> s]
-  real, dimension(:,:),             optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage
+  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)       :: frac_shelf_h !< Fractional ice shelf coverage [nondim]
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
@@ -342,9 +342,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
   nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
 
   ice_shelf = .false.
-  if (present(frac_shelf_h)) then
-    if (associated(frac_shelf_h)) ice_shelf = .true.
-  endif
+  if (present(frac_shelf_h)) ice_shelf = .true.
 
   if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90")
 
@@ -621,7 +619,7 @@ subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h
   real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h     !< Current 3D grid obtained after the
                                                                      !! last time step [H ~> m or kg-2]
   logical,                       optional, intent(in)    :: debug    !< If true, show the call tree
-  real, dimension(:,:),          optional, pointer       :: frac_shelf_h !< Fractional ice shelf coverage
+  real, dimension(SZI_(G),SZJ_(G)),  optional, intent(in):: frac_shelf_h !< Fractional ice shelf coverage [nondim]
   ! Local variables
   integer :: nk, i, j, k
   real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
@@ -633,7 +631,7 @@ subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h
   if (show_call_tree) call callTree_enter("ALE_build_grid(), MOM_ALE.F90")
   use_ice_shelf = .false.
   if (present(frac_shelf_h)) then
-    if (associated(frac_shelf_h)) use_ice_shelf = .true.
+    use_ice_shelf = .true.
   endif
 
   ! Build new grid. The new grid is stored in h_new. The old grid is h.