Merge branch 'origin/offline_tracers' of https://github.com/ashao/MOM6 …

…into ashao-origin/offline_tracers

config_src/solo_driver/MOM_driver.F90

@@ -47,6 +47,7 @@ program MOM_main
   use MOM_diag_mediator,   only : diag_mediator_close_registration, diag_mediator_end
   use MOM,                 only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
   use MOM,                 only : calculate_surface_state, finish_MOM_initialization
+  use MOM,                 only : step_tracers
   use MOM_domains,         only : MOM_infra_init, MOM_infra_end
   use MOM_error_handler,   only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe
   use MOM_error_handler,   only : callTree_enter, callTree_leave, callTree_waypoint
@@ -134,6 +135,7 @@ program MOM_main
                                   ! representation of time_step.
   real :: time_step               ! The time step of a call to step_MOM in seconds.
   real :: dt                      ! The baroclinic dynamics time step, in seconds.
+  real :: dt_off                  ! Offline time step in seconds
   integer :: ntstep               ! The number of baroclinic dynamics time steps
                                   ! within time_step.
 
@@ -169,6 +171,14 @@ program MOM_main
   logical :: unit_in_use
   integer :: initClock, mainClock, termClock
 
+  logical :: do_online          ! If true, use the model in prognostic mode where
+                                ! the barotropic and baroclinic dynamics, thermodynamics,
+                                ! etc. are stepped forward integrated in time.
+                                ! If false, then all of the above are bypassed with all
+                                ! fields necessary to integrate only the tracer advection
+                                ! and diffusion equation are read in from files stored from
+                                ! a previous integration of the prognostic model
+
   type(MOM_control_struct),  pointer :: MOM_CSp => NULL()
   type(surface_forcing_CS),  pointer :: surface_forcing_CSp => NULL()
   type(sum_output_CS),       pointer :: sum_output_CSp => NULL()
@@ -250,17 +260,18 @@ program MOM_main
   else
     Start_time = set_time(0,days=0)
   endif
-
+  
   if (sum(date) >= 0) then
     ! 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, param_file, dirs, MOM_CSp, segment_start_time)
+    ! Note the not before CS%d
+    call initialize_MOM(Time, param_file, dirs, MOM_CSp, segment_start_time, do_online_out = do_online)
   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
-    call initialize_MOM(Time, param_file, dirs, MOM_CSp)
+    call initialize_MOM(Time, param_file, dirs, MOM_CSp, do_online_out=do_online)
   endif
   fluxes%C_p = MOM_CSp%tv%C_p  ! Copy the heat capacity for consistency.
 
@@ -298,8 +309,12 @@ program MOM_main
   call get_param(param_file, mod, "DT_FORCING", time_step, &
                  "The time step for changing forcing, coupling with other \n"//&
                  "components, or potentially writing certain diagnostics. \n"//&
-                 "The default value is given by DT.", units="s", default=dt)
-
+                 "The default value is given by DT.", units="s", default=dt)             
+  if (.not. do_online) then
+    call get_param(param_file, mod, "DT_OFFLINE", time_step, &
+                   "Time step for the offline time step")
+    dt = time_step
+  endif
   ntstep = MAX(1,ceiling(time_step/dt - 0.001))
 
   Time_step_ocean = set_time(int(floor(time_step+0.5)))
@@ -350,6 +365,7 @@ program MOM_main
                  "The interval in units of TIMEUNIT between saves of the \n"//&
                  "energies of the run and other globally summed diagnostics.", &
                  default=set_time(int(time_step+0.5)), timeunit=Time_unit)
+
   call log_param(param_file, mod, "ELAPSED TIME AS MASTER", elapsed_time_master)
 
   ! Close the param_file.  No further parsing of input is possible after this.
@@ -399,8 +415,10 @@ program MOM_main
     call callTree_enter("Main loop, MOM_driver.F90",n)
 
     ! Set the forcing for the next steps.
-    call set_forcing(state, fluxes, Time, Time_step_ocean, grid, &
+    if (do_online) then
+        call set_forcing(state, fluxes, Time, Time_step_ocean, grid, &
                      surface_forcing_CSp)
+    endif
     if (MOM_CSp%debug) then
       call MOM_forcing_chksum("After set forcing", fluxes, grid, haloshift=0)
     endif
@@ -423,7 +441,8 @@ program MOM_main
 
     ! This call steps the model over a time time_step.
     Time1 = Master_Time ; Time = Master_Time
-    call step_MOM(fluxes, state, Time1, time_step, MOM_CSp)
+    if (do_online)        call step_MOM(fluxes, state, Time1, time_step, MOM_CSp)
+    if (.not. do_online)  call step_tracers(fluxes, state, Time1, time_step, MOM_CSp)
 
 !    Time = Time + Time_step_ocean
 !  This is here to enable fractional-second time steps.
@@ -450,17 +469,21 @@ program MOM_main
                             surface_forcing_CSp%handles)
     call disable_averaging(MOM_CSp%diag)
 
-    if (fluxes%fluxes_used) then
-      call enable_averaging(fluxes%dt_buoy_accum, Time, MOM_CSp%diag)
-      call forcing_diagnostics(fluxes, state, fluxes%dt_buoy_accum, grid, &
-                               MOM_CSp%diag, surface_forcing_CSp%handles)
-      call accumulate_net_input(fluxes, state, fluxes%dt_buoy_accum, grid, sum_output_CSp)
-      call disable_averaging(MOM_CSp%diag)
-    else
-      call MOM_error(FATAL, "The solo MOM_driver is not yet set up to handle "//&
-             "thermodynamic time steps that are longer than the couping timestep.")
+    if (do_online) then
+      if (fluxes%fluxes_used) then
+        call enable_averaging(fluxes%dt_buoy_accum, Time, MOM_CSp%diag)
+        call forcing_diagnostics(fluxes, state, fluxes%dt_buoy_accum, grid, &
+                                 MOM_CSp%diag, surface_forcing_CSp%handles)
+        call accumulate_net_input(fluxes, state, fluxes%dt_buoy_accum, grid, sum_output_CSp)
+        call disable_averaging(MOM_CSp%diag)
+      else
+        call MOM_error(FATAL, "The solo MOM_driver is not yet set up to handle "//&
+               "thermodynamic time steps that are longer than the coupling timestep.")
+      endif
     endif
 
+
+
 !  See if it is time to write out the energy.
     if ((Time + (Time_step_ocean/2) > write_energy_time) .and. &
         (MOM_CSp%dt_trans == 0.0)) then
@@ -502,7 +525,7 @@ program MOM_main
   if (Restart_control>=0) then
     if (MOM_CSp%dt_trans > 0.0) call MOM_error(WARNING, "End of MOM_main reached "//&
          "with a non-zero dt_trans.  Additional restart fields are required.")
-     if (.not.fluxes%fluxes_used) call MOM_error(FATAL, "End of MOM_main reached "//&
+     if (.not.fluxes%fluxes_used .and. do_online) call MOM_error(FATAL, "End of MOM_main reached "//&
          "with unused buoyancy fluxes.  For conservation, the ocean restart "//&
          "files can only be created after the buoyancy forcing is applied.")
 

src/ALE/MOM_ALE.F90

@@ -97,14 +97,17 @@ module MOM_ALE
   integer, dimension(:), allocatable :: id_tracer_remap_tendency      !< diagnostic id 
   integer, dimension(:), allocatable :: id_Htracer_remap_tendency     !< diagnostic id 
   integer, dimension(:), allocatable :: id_Htracer_remap_tendency_2d  !< diagnostic id 
-  logical, dimension(:), allocatable :: do_tendency_diag              !< flag for doing diagnostics 
+  logical, dimension(:), allocatable :: do_tendency_diag              !< flag for doing diagnostics
+  integer                            :: id_dzRegrid
 
 end type
 
 ! Publicly available functions
 public ALE_init
 public ALE_end
-public ALE_main 
+public ALE_main
+public ALE_main_offline
+public ALE_offline_tracer_final
 public ALE_build_grid
 public ALE_remap_scalar
 public pressure_gradient_plm
@@ -282,6 +285,9 @@ subroutine ALE_register_diags(Time, G, diag, C_p, Reg, CS)
   CS%id_Htracer_remap_tendency(:)    = -1
   CS%id_Htracer_remap_tendency_2d(:) = -1
 
+  CS%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,Time, &
+      'Change in interface height due to ALE regridding', 'meter')
+
   if(ntr > 0) then 
 
     do m=1,ntr
@@ -375,7 +381,6 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt)
   type(tracer_registry_type),                 pointer       :: Reg !< Tracer registry structure
   type(ALE_CS),                               pointer       :: CS  !< Regridding parameters and options
   real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main()
-
   ! 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)) :: h_new ! New 3D grid obtained after last time step (m or Pa)
@@ -399,6 +404,7 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt)
   if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)")
 
   ! Remap all variables from old grid h onto new grid h_new
+
   call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, -dzRegrid, Reg, &
                              u, v, CS%show_call_tree, dt )
 
@@ -414,8 +420,114 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt)
   enddo
 
   if (CS%show_call_tree) call callTree_leave("ALE_main()")
+
+  if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag)
+
+
 end subroutine ALE_main
 
+!> Takes care of (1) building a new grid and (2) remapping all variables between
+!! the old grid and the new grid. The creation of the new grid can be based
+!! on z coordinates, target interface densities, sigma coordinates or any
+!! arbitrary coordinate system.
+subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, dt)
+  type(ocean_grid_type),                      intent(in)    :: G   !< Ocean grid informations
+  type(verticalGrid_type),                    intent(in)    :: GV  !< Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: h   !< Current 3D grid obtained after last time step (m or Pa)
+  type(thermo_var_ptrs),                      intent(inout) :: tv  !< Thermodynamic variable structure
+  type(tracer_registry_type),                 pointer       :: Reg !< Tracer registry structure
+  type(ALE_CS),                               pointer       :: CS  !< Regridding parameters and options
+  real,                             optional, intent(in)    :: dt  !< Time step between calls to ALE_main()
+  ! 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)) :: h_new ! New 3D grid obtained after last time step (m or Pa)
+  integer :: nk, i, j, k, isc, iec, jsc, jec
+
+  nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
+
+  if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90")
+
+  if (present(dt)) then
+    call ALE_update_regrid_weights( dt, CS )
+  endif
+  dzRegrid(:,:,:) = 0.0
+
+  ! Build new grid. The new grid is stored in h_new. The old grid is h.
+  ! Both are needed for the subsequent remapping of variables.
+  call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid )
+
+  call check_grid( G, GV, h, 0. )
+
+  if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)")
+
+  ! Remap all variables from old grid h onto new grid h_new
+
+  call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, -dzRegrid, Reg, &
+                             debug=CS%show_call_tree, dt=dt )
+
+  if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)")
+
+  ! Override old grid with new one. The new grid 'h_new' is built in
+  ! one of the 'build_...' routines above.
+!$OMP parallel do default(none) shared(isc,iec,jsc,jec,nk,h,h_new,CS)
+  do k = 1,nk
+    do j = jsc-1,jec+1 ; do i = isc-1,iec+1
+      h(i,j,k) = h_new(i,j,k)
+    enddo ; enddo
+  enddo
+
+  if (CS%show_call_tree) call callTree_leave("ALE_main()")
+  if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag)
+
+end subroutine ALE_main_offline
+
+!> Remaps all tracers from h onto h_target. This is intended to be called when tracers
+!! are done offline. In the case where transports don't quite conserve, we still want to
+!! make sure that layer thicknesses offline do not drift too far away from the online model
+subroutine ALE_offline_tracer_final( G, GV, h, h_target, Reg, CS)
+  type(ocean_grid_type),                      intent(in)    :: G   !< Ocean grid informations
+  type(verticalGrid_type),                    intent(in)    :: GV  !< Ocean vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: h   !< Current 3D grid obtained after last time step (m or Pa)
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(in)    :: h_target   !< Current 3D grid obtained after last time step (m or Pa)
+  type(tracer_registry_type),                 pointer       :: Reg !< Tracer registry structure
+  type(ALE_CS),                               pointer       :: CS  !< Regridding parameters and options
+  ! Local variables
+
+  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
+  integer :: nk, i, j, k, isc, iec, jsc, jec
+
+  nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
+
+  if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90")
+
+  ! It does not seem that remap_all_state_vars uses dzRegrid for tracers, only for u, v
+  dzRegrid(:,:,:) = 0.0
+
+  call check_grid( G, GV, h, 0. )
+  call check_grid( G, GV, h_target, 0. )
+
+  if (CS%show_call_tree) call callTree_waypoint("Source and target grids checked (ALE_offline_tracer)")
+
+  ! Remap all variables from old grid h onto new grid h_new
+
+  call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_target, -dzRegrid, Reg, &
+                             debug=CS%show_call_tree )
+
+  if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_offline_tracer)")
+
+  ! Override old grid with new one. The new grid 'h_new' is built in
+  ! one of the 'build_...' routines above.
+!$OMP parallel do default(none) shared(isc,iec,jsc,jec,nk,h,h_target,CS)
+  do k = 1,nk
+    do j = jsc-1,jec+1 ; do i = isc-1,iec+1
+      h(i,j,k) = h_target(i,j,k)
+    enddo ; enddo
+  enddo
+
+  if (CS%show_call_tree) call callTree_leave("ALE_offline_tracer()")
+
+end subroutine ALE_offline_tracer_final
+
 !> Check grid for negative thicknesses
 subroutine check_grid( G, GV, h, threshold )
   type(ocean_grid_type),                     intent(in) :: G !< Ocean grid structure