Merge branch 'nicjhan-62-z-remapping-uv-grids' into dev/master

config_src/coupled_driver/ocean_model_MOM.F90

@@ -37,7 +37,7 @@ module ocean_model_mod
 use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
 use MOM, only : calculate_surface_state
 use MOM_constants, only : CELSIUS_KELVIN_OFFSET
-use MOM_diag_mediator, only : enable_averaging, disable_averaging
+use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging
 use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end
 use MOM_domains, only : pass_vector, AGRID, BGRID_NE, CGRID_NE
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
@@ -502,6 +502,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time)
   type(ocean_public_type),           intent(inout) :: Ocean_sfc
   type(ocean_state_type),            pointer       :: Ocean_state
   type(time_type),                   intent(in)    :: Time
+
 !   This subroutine terminates the model run, saving the ocean state in a
 ! restart file and deallocating any data associated with the ocean.
 
@@ -512,7 +513,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time)
 !  (in)      Time - The model time, used for writing restarts.
 
   call ocean_model_save_restart(Ocean_state, Time)
-  call diag_mediator_end(Time)
+  call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag)
   call MOM_end(Ocean_state%MOM_CSp)
   if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp)
 end subroutine ocean_model_end

config_src/ice_solo_driver/ice_shelf_driver.F90

@@ -416,7 +416,7 @@ program SHELF_main
     close(unit)
   endif
 
-  call diag_mediator_end(Time, end_diag_manager=.true.)
+  call diag_mediator_end(Time, ice_shelf_CSp%diag, end_diag_manager=.true.)
   call cpu_clock_end(termClock)
 
   call io_infra_end ; call MOM_infra_end

config_src/solo_driver/MOM_driver.F90

@@ -517,7 +517,7 @@ program MOM_main
   endif
 
   call callTree_waypoint("End MOM_main")
-  call diag_mediator_end(Time, end_diag_manager=.true.)
+  call diag_mediator_end(Time, MOM_CSp%diag, end_diag_manager=.true.)
   call cpu_clock_end(termClock)
 
   call io_infra_end ; call MOM_infra_end

src/ALE/MOM_remapping.F90

@@ -534,15 +534,22 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 )
     if (k<=n0) then; hTmp = h0(k); else; hTmp = 0.; endif
     z0 = z0 + hTmp ; z1 = z1 + ( hTmp + ( dx(k+1) - dx(k) ) )
   enddo
+  ! Maximum error based on guess at maximum roundoff
   if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then
-    write(0,*) 'h0=',h0
-    write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1)
-    write(0,*) 'u0=',u0
-    write(0,*) 'u1=',u1
-    write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0
-    write(0,*) 'err0,errF=',err0,err2
-    call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//&
-         'Total stuff on h0 and hF differ by more than roundoff' )
+    ! Maximum relative error
+    if (abs(totalHU2-totalHU0) / totalHU2 > 1e-09) then
+      ! Maximum absolute error
+      if (abs(totalHU2-totalHU0) > 1e-18) then
+        write(0,*) 'h0=',h0
+        write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1)
+        write(0,*) 'u0=',u0
+        write(0,*) 'u1=',u1
+        write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0
+        write(0,*) 'err0,errF=',err0,err2
+        call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//&
+             'Total stuff on h0 and hF differ by more than maximum errors' )
+      endif
+    endif
   endif
 #endif
 

src/core/MOM.F90

@@ -337,7 +337,8 @@ module MOM
 use MOM_cpu_clock,            only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
 use MOM_cpu_clock,            only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
 use MOM_coms,                 only : reproducing_sum
-use MOM_diag_mediator,        only : diag_mediator_init, enable_averaging, diag_set_thickness_ptr
+use MOM_diag_mediator,        only : diag_mediator_init, enable_averaging
+use MOM_diag_mediator,        only : diag_set_thickness_ptr, diag_update_target_grids
 use MOM_diag_mediator,        only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator,        only : register_diag_field, register_static_field
 use MOM_diag_mediator,        only : register_scalar_field
@@ -977,6 +978,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         call do_group_pass(CS%pass_uv_T_S_h, G%Domain)
         call cpu_clock_end(id_clock_pass)
 
+        ! Whenever thickness changes let the diag manager know, target grids
+        ! for vertical remapping may need to be regenerated.
+        call diag_update_target_grids(CS%diag)
+
         if (CS%debug) then
           call uchksum(u,"Post-dia first u", G, haloshift=2)
           call vchksum(v,"Post-dia first v", G, haloshift=2)
@@ -1039,6 +1044,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         else
           dtth = dt*min(ntstep,n_max-n+1)
         endif
+
         call enable_averaging(dtth,Time_local+set_time(int(floor(dtth-dt+0.5))), CS%diag)
         call cpu_clock_begin(id_clock_thick_diff)
         if (associated(CS%VarMix)) call calc_slope_functions(h, CS%tv, dt, G, CS%VarMix)
@@ -1050,6 +1056,11 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
         call cpu_clock_end(id_clock_pass)
         call disable_averaging(CS%diag)
         if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)")
+
+        ! Whenever thickness changes let the diag manager know, target grids
+        ! for vertical remapping may need to be regenerated.
+        call diag_update_target_grids(CS%diag)
+
       endif
     endif
 
@@ -1159,11 +1170,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
     if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, &
                                            CS%visc, dt, G, CS%MEKE_CSp)
-
     call disable_averaging(CS%diag)
     call cpu_clock_end(id_clock_dynamics)
 
-
     CS%dt_trans = CS%dt_trans + dt
     if (thermo_does_span_coupling) then
       do_advection = (CS%dt_trans + 0.5*dt > dt_therm)
@@ -1270,7 +1279,12 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
             call hchksum(CS%tv%S,"Post-ALE S", G, haloshift=1)
             call check_redundant("Post-ALE ", u, v, G)
           endif
-        endif   
+        endif
+
+        ! Whenever thickness changes let the diag manager know, target grids
+        ! for vertical remapping may need to be regenerated. This needs to
+        ! happen after the H update and before the next post_data.
+        call diag_update_target_grids(CS%diag)
 
         call cpu_clock_begin(id_clock_pass)
         call do_group_pass(CS%pass_uv_T_S_h, G%Domain)
@@ -1928,10 +1942,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   endif
   call callTree_waypoint("state variables allocated (initialize_MOM)")
 
-  ! Set up a pointers h within diag mediator control structure,
-  ! this needs to occur _after_ CS%h has been allocated.
-  call diag_set_thickness_ptr(CS%h, diag)
-
   ! Set the fields that are needed for bitwise identical restarting
   ! the time stepping scheme.
   call restart_init(G, param_file, CS%restart_CSp)
@@ -1980,11 +1990,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   call cpu_clock_end(id_clock_MOM_init)
   call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")
 
-  ! Initialize the diagnostics mask arrays. 
-  ! This step has to be done after call to MOM_initialize_state
-  ! and before MOM_diagnostics_init
-  call diag_masks_set(G, CS%missing, diag)
-
   if (CS%use_ALE_algorithm) then
     ! For now, this has to follow immediately after MOM_initialize_state because
     ! the call to initialize_ALE can change CS%h, etc.  initialize_ALE should
@@ -2008,9 +2013,24 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
     endif
   endif
 
-  ! This call sets up the diagnostic axes.
-  call cpu_clock_begin(id_clock_MOM_init)
+  ! Initialize the diagnostics mask arrays.
+  ! This step has to be done after call to MOM_initialize_state
+  ! and before MOM_diagnostics_init
+  call diag_masks_set(G, CS%missing, diag)
+
+  ! Set up a pointers h within diag mediator control structure,
+  ! this needs to occur _after_ CS%h has been allocated.
+  call diag_set_thickness_ptr(CS%h, diag)
+
+  ! This call sets up the diagnostic axes. These are needed,
+  ! e.g. to generate the target grids below.
   call set_axes_info(G, param_file, diag)
+
+  ! Whenever thickness changes let the diag manager know, target grids
+  ! for vertical remapping may need to be regenerated. This needs to
+  call diag_update_target_grids(diag)
+
+  call cpu_clock_begin(id_clock_MOM_init)
   if (CS%use_ALE_algorithm) then
     call ALE_writeCoordinateFile( CS%ALE_CSp, G, dirs%output_directory )
   endif

src/core/MOM_dynamics_legacy_split.F90

@@ -78,7 +78,7 @@ module MOM_dynamics_legacy_split
 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator, only : register_diag_field, register_static_field
-use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl
+use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids
 use MOM_domains, only : MOM_domains_init, pass_var, pass_vector
 use MOM_domains, only : pass_var_start, pass_var_complete
 use MOM_domains, only : pass_vector_start, pass_vector_complete
@@ -631,8 +631,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
     if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then
       call cpu_clock_begin(id_clock_continuity)
       call continuity(u, v, h, hp, uh_in, vh_in, dt, G, &
-                      CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, &
-                      visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
+                      CS%continuity_CSp, OBC=CS%OBC, &
+                      visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, &
+                      BT_cont=CS%BT_cont)
       call cpu_clock_end(id_clock_continuity)
       if (BT_cont_BT_thick) then
         call cpu_clock_begin(id_clock_pass)
@@ -716,8 +717,8 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
 ! hp = h + dt * div . uh
   call cpu_clock_begin(id_clock_continuity)
   call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, &
-                  CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, &
-                  u_av, v_av, BT_cont=CS%BT_cont)
+                  CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, &
+                  CS%visc_rem_v, u_av, v_av, BT_cont=CS%BT_cont)
   call cpu_clock_end(id_clock_continuity)
 
   call cpu_clock_begin(id_clock_pass)
@@ -931,6 +932,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
                     CS%visc_rem_u, CS%visc_rem_v, u_av, v_av, &
                     uhbt_out, vhbt_out, u, v)
     call cpu_clock_end(id_clock_continuity)
+    ! Whenever thickness changes let the diag manager know, target grids
+    ! for vertical remapping may need to be regenerated.
+    call diag_update_target_grids(CS%diag)
     if (G%nonblocking_updates) then
       call cpu_clock_begin(id_clock_pass)
       pid_h = pass_var_start(h, G%Domain)
@@ -969,6 +973,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, &
                     CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, &
                     CS%visc_rem_u, CS%visc_rem_v, u_av, v_av)
     call cpu_clock_end(id_clock_continuity)
+    ! Whenever thickness changes let the diag manager know, target grids
+    ! for vertical remapping may need to be regenerated.
+    call diag_update_target_grids(CS%diag)
     call cpu_clock_begin(id_clock_pass)
     call pass_var(h, G%Domain)
     call cpu_clock_end(id_clock_pass)

src/core/MOM_dynamics_split_RK2.F90

@@ -78,7 +78,7 @@ module MOM_dynamics_split_RK2
 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator, only : register_diag_field, register_static_field
-use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl
+use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids
 use MOM_domains, only : MOM_domains_init
 use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR
 use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
@@ -914,6 +914,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
                   CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, &
                   CS%visc_rem_u, CS%visc_rem_v, u_av, v_av)
   call cpu_clock_end(id_clock_continuity)
+  ! Whenever thickness changes let the diag manager know, target grids
+  ! for vertical remapping may need to be regenerated.
+  call diag_update_target_grids(CS%diag)
   call cpu_clock_begin(id_clock_pass)
   call do_group_pass(CS%pass_h, G%Domain)
   call cpu_clock_end(id_clock_pass)

src/core/MOM_dynamics_unsplit.F90

@@ -78,7 +78,7 @@ module MOM_dynamics_unsplit
 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging
 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator, only : register_diag_field, register_static_field
-use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl
+use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids
 use MOM_domains, only : MOM_domains_init, pass_var, pass_vector
 use MOM_domains, only : pass_var_start, pass_var_complete
 use MOM_domains, only : pass_vector_start, pass_vector_complete
@@ -268,7 +268,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
 ! uh = u*h
 ! hp = h + dt/2 div . uh
   call cpu_clock_begin(id_clock_continuity)
-  call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, OBC=CS%OBC)
+  call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, &
+                  OBC=CS%OBC)
   call cpu_clock_end(id_clock_continuity)
   call cpu_clock_begin(id_clock_pass)
   call pass_var(hp, G%Domain)
@@ -438,6 +439,9 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
   call continuity(upp, vpp, hp, h, uh, vh, &
                   (dt*0.5), G, CS%continuity_CSp, OBC=CS%OBC)
   call cpu_clock_end(id_clock_continuity)
+  ! Whenever thickness changes let the diag manager know, target grids
+  ! for vertical remapping may need to be regenerated.
+  call diag_update_target_grids(CS%diag)
   call cpu_clock_begin(id_clock_pass)
   call pass_var(h, G%Domain)
   call pass_vector(uh, vh, G%Domain)

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -282,7 +282,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   call cpu_clock_begin(id_clock_continuity)
   ! This is a duplicate caclulation of the last continuity from the previous step
   ! and could/should be optimized out. -AJA
-  call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, OBC=CS%OBC)
+  call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, &
+                  OBC=CS%OBC)
   call cpu_clock_end(id_clock_continuity)
   call cpu_clock_begin(id_clock_pass)
   call pass_var(hp, G%Domain)