Add option to use QG Leith with biharmonic

src/ALE/MOM_regridding.F90

@@ -799,7 +799,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_
       call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h )
       call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
     case ( REGRIDDING_SIGMA_SHELF_ZSTAR)
-      call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h )
+      call build_zstar_grid( CS, G, GV, h, dzInterface )
       call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
     case ( REGRIDDING_SIGMA )
       call build_sigma_grid( CS, G, GV, h, dzInterface )

src/core/MOM.F90

@@ -139,6 +139,8 @@ module MOM
 ! ODA modules
 use MOM_oda_driver_mod,        only : ODA_CS, oda, init_oda, oda_end
 use MOM_oda_driver_mod,        only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments
+use MOM_oda_incupd,            only : oda_incupd_CS, init_oda_incupd_diags
+
 ! Offline modules
 use MOM_offline_main,          only : offline_transport_CS, offline_transport_init, update_offline_fields
 use MOM_offline_main,          only : insert_offline_main, extract_offline_main, post_offline_convergence_diags
@@ -369,6 +371,8 @@ module MOM
   type(sponge_CS),               pointer :: sponge_CSp => NULL()
     !< Pointer to the layered-mode sponge control structure
   type(ALE_sponge_CS),           pointer :: ALE_sponge_CSp => NULL()
+    !< Pointer to the oda incremental update control structure
+  type(oda_incupd_CS),           pointer :: oda_incupd_CSp => NULL()
     !< Pointer to the ALE-mode sponge control structure
   type(ALE_CS),                  pointer :: ALE_CSp => NULL()
     !< Pointer to the Arbitrary Lagrangian Eulerian (ALE) vertical coordinate control structure
@@ -1667,6 +1671,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   type(ocean_OBC_type), pointer :: OBC_in => NULL()
   type(sponge_CS), pointer :: sponge_in_CSp => NULL()
   type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => NULL()
+  type(oda_incupd_CS),pointer :: oda_incupd_in_CSp => NULL()
 
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
@@ -2425,12 +2430,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
       call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in)
       call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
-          sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in, &
+          sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, &
           frac_shelf_h=frac_shelf_in)
     else
       call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
-          sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in)
+          sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in)
     endif
 
     if (use_temperature) then
@@ -2472,12 +2477,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
       call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h)
       call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
-          CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in, &
+          CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, &
           frac_shelf_h=CS%frac_shelf_h)
     else
       call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
-          CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in)
+          CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in)
     endif
   endif
 
@@ -2661,7 +2666,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   else
     call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, &
                               CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, &
-                              CS%sponge_CSp, CS%ALE_sponge_CSp)
+                              CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp)
   endif
 
   if (associated(CS%sponge_CSp)) &
@@ -2670,6 +2675,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   if (associated(CS%ALE_sponge_CSp)) &
     call init_ALE_sponge_diags(Time, G, diag, CS%ALE_sponge_CSp, US)
 
+  if (associated(CS%oda_incupd_CSp)) &
+    call init_oda_incupd_diags(Time, G, GV, diag, CS%oda_incupd_CSp, US)
+
+
   call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp)
   call tracer_hor_diff_init(Time, G, GV, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, &
                             CS%tracer_diff_CSp)

src/framework/MOM_diag_mediator.F90

@@ -1969,10 +1969,10 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
                                                          !! integrated). Default/absent for intensive.
   ! Local variables
   real :: MOM_missing_value
-  type(diag_ctrl), pointer :: diag_cs => NULL()
-  type(axes_grp), pointer :: remap_axes => null()
-  type(axes_grp), pointer :: axes => null()
-  type(axes_grp), pointer :: axes_d2 => null()
+  type(diag_ctrl), pointer :: diag_cs
+  type(axes_grp), pointer :: remap_axes
+  type(axes_grp), pointer :: axes
+  type(axes_grp), pointer :: axes_d2
   integer :: dm_id, i, dl
   character(len=256) :: msg, cm_string
   character(len=256) :: new_module_name
@@ -2097,8 +2097,8 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
 
     new_module_name = trim(module_name)//'_d2'
 
+    axes_d2 => null()
     if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then
-      axes_d2 => null()
       if (axes_in%id == diag_cs%axesTL%id) then
         axes_d2 => diag_cs%dsamp(dl)%axesTL
       elseif (axes_in%id == diag_cs%axesBL%id) then
@@ -2129,6 +2129,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
               //trim(new_module_name)//"-"//trim(field_name))
       endif
     endif
+
     ! Register the native diagnostic
     if (associated(axes_d2)) then
        active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, &

src/initialization/MOM_state_initialization.F90

@@ -11,7 +11,7 @@ module MOM_state_initialization
 use MOM_cpu_clock, only :  CLOCK_ROUTINE, CLOCK_LOOP
 use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
 use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID
-use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
+use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, NOTE, WARNING, is_root_pe
 use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser, only : get_param, read_param, log_param, param_file_type
 use MOM_file_parser, only : log_version
@@ -94,6 +94,9 @@ module MOM_state_initialization
 use MOM_remapping, only : remapping_CS, initialize_remapping
 use MOM_remapping, only : remapping_core_h
 use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer
+use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd
+use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field
+use MOM_oda_incupd, only: calc_oda_increments, output_oda_incupd_inc
 
 implicit none ; private
 
@@ -114,7 +117,7 @@ module MOM_state_initialization
 !! conditions or by reading them from a restart (or saves) file.
 subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                 restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
-                                ALE_sponge_CSp, OBC, Time_in, frac_shelf_h)
+                                ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h)
   type(ocean_grid_type),      intent(inout) :: G    !< The ocean's grid structure.
   type(verticalGrid_type),    intent(in)    :: GV   !< The ocean's vertical grid structure.
   type(unit_scale_type),      intent(in)    :: US   !< A dimensional unit scaling type
@@ -140,6 +143,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   type(sponge_CS),            pointer       :: sponge_CSp !< The layerwise sponge control structure.
   type(ALE_sponge_CS),        pointer       :: ALE_sponge_CSp !< The ALE sponge control structure.
   type(ocean_OBC_type),       pointer       :: OBC   !< The open boundary condition control structure.
+  type(oda_incupd_CS),        pointer       :: oda_incupd_CSp !< The oda_incupd control structure.
   type(time_type), optional,  intent(in)    :: Time_in !< Time at the start of the run segment.
   real, dimension(SZI_(G),SZJ_(G)), &
                      optional, intent(in)   :: frac_shelf_h    !< The fraction of the grid cell covered
@@ -157,7 +161,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   logical :: from_Z_file, useALE
   logical :: new_sim
   integer :: write_geom
-  logical :: use_temperature, use_sponge, use_OBC
+  logical :: use_temperature, use_sponge, use_OBC, use_oda_incupd
   logical :: use_EOS     ! If true, density is calculated from T & S using an equation of state.
   logical :: depress_sfc ! If true, remove the mass that would be displaced
                          ! by a large surface pressure by squeezing the column.
@@ -478,6 +482,16 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                   dt=dt, initial=.true.)
     endif
   endif
+
+  ! Initialized assimilative incremental update (oda_incupd) structure and
+  ! register restart.
+  call get_param(PF, mdl, "ODA_INCUPD", use_oda_incupd, &
+                 "If true, oda incremental updates will be applied "//&
+                 "everywhere in the domain.", default=.false.)
+  if (use_oda_incupd) then
+    call initialize_oda_incupd_fixed(G, GV, US, oda_incupd_CSp, restart_CS)
+  endif
+
   ! This is the end of the block of code that might have initialized fields
   ! internally at the start of a new run.
 
@@ -615,6 +629,13 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h)
   call callTree_leave('MOM_initialize_state()')
 
+  ! Set-up of data Assimilation with incremental update
+  if (use_oda_incupd) then
+    call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, &
+                                    PF, oda_incupd_CSp, restart_CS, Time)
+  endif
+
+
 end subroutine MOM_initialize_state
 
 !> Reads the layer thicknesses or interface heights from a file.
@@ -2022,6 +2043,157 @@ end function separate_idamp_for_uv
 
 end subroutine initialize_sponges_file
 
+subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, &
+                                      oda_incupd_CSp, restart_CS, Time)
+  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
+  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
+  type(unit_scale_type),   intent(in)    :: US  !< A dimensional unit scaling type
+  logical,                 intent(in)    :: use_temperature !< If true, T & S are state variables.
+  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic
+                                              !! variables.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                             intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2] (in)
+
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
+                             intent(in) :: u    !< The zonal velocity that is being
+                                                 !! initialized [L T-1 ~> m s-1]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
+                             intent(in) :: v    !< The meridional velocity that is being
+                                                !! initialized [L T-1 ~> m s-1]
+  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
+  type(oda_incupd_CS),     pointer    :: oda_incupd_CSp  !< A pointer that is set to point to the control
+                                                  !! structure for this module.
+  type(MOM_restart_CS),    pointer    :: restart_CS !< A pointer to the restart control
+                                                    !! structure.
+  type(time_type),         intent(in) :: Time !< Time at the start of the run segment. Time_in
+                                              !! overrides any value set for
+                                              !Time.
+  ! Local variables
+  real, allocatable, dimension(:,:,:) :: hoda ! The layer thk inc. and oda layer thk [H ~> m or kg m-2].
+  real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields
+  real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading oda fields
+
+  integer :: i, j, k, is, ie, js, je, nz
+  integer :: isd, ied, jsd, jed
+
+  integer, dimension(4) :: siz
+  integer :: nz_data  ! The size of the sponge source grid
+  logical :: oda_inc  ! input files are increments (true) or full fields (false)
+  logical :: save_inc ! save increments if using full fields
+  logical :: uv_inc   ! use u and v increments
+  logical :: reset_ncount ! reset ncount to zero if true
+
+  character(len=40)  :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var
+  character(len=40)  :: mdl = "initialize_oda_incupd_file"
+  character(len=200) :: inc_file, uv_inc_file  ! Strings for filenames
+  character(len=200) :: filename, inputdir ! Strings for file/path and path.
+
+!  logical :: use_ALE ! True if ALE is being used, False if in layered mode
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+
+  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
+  inputdir = slasher(inputdir)
+
+  call get_param(param_file, mdl, "ODA_INCUPD_FILE", inc_file, &
+                 "The name of the file with the T,S,h increments.", &
+                 fail_if_missing=.true.)
+  call get_param(param_file, mdl, "ODA_INCUPD_INC", oda_inc, &
+                 "INCUPD files are increments and not full fields.", &
+                 default=.true.)
+  if (.not.oda_inc) then
+    call get_param(param_file, mdl, "ODA_INCUPD_SAVE", save_inc, &
+                   "If true, save the increments when using full fields.", &
+                   default=.false.)
+  endif
+  call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, &
+                 "If True, reinitialize number of updates already done, ncount.",&
+                 default=.true.)
+  if (.not.oda_inc .and. .not.reset_ncount) &
+    call MOM_error(FATAL, " initialize_oda_incupd: restarting during update "// &
+                 "necessitates increments input file")
+
+  call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, &
+                 "The name of the potential temperature inc. variable in "//&
+                 "ODA_INCUPD_FILE.", default="ptemp_inc")
+  call get_param(param_file, mdl, "ODA_SALTINC_VAR", salinc_var, &
+                 "The name of the salinity inc. variable in "//&
+                 "ODA_INCUPD_FILE.", default="sal_inc")
+  call get_param(param_file, mdl, "ODA_THK_VAR", h_var, &
+                 "The name of the layer thickness variable in "//&
+                 "ODA_INCUPD_FILE.", default="h")
+  call get_param(param_file, mdl, "ODA_INCUPD_UV", uv_inc, &
+                 "use U,V increments.", &
+                 default=.true.)
+  call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, &
+                 "The name of the file with the U,V increments.", &
+                 default=inc_file)
+  call get_param(param_file, mdl, "ODA_UINC_VAR", uinc_var, &
+                 "The name of the zonal vel. inc. variable in "//&
+                 "ODA_INCUPD_FILE.", default="u_inc")
+  call get_param(param_file, mdl, "ODA_VINC_VAR", vinc_var, &
+                 "The name of the meridional vel. inc. variable in "//&
+                 "ODA_INCUPD_FILE.", default="v_inc")
+
+!  call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.)
+
+  ! Read in incremental update for tracers
+  filename = trim(inputdir)//trim(inc_file)
+  call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_FILE", filename)
+  if (.not.file_exists(filename, G%Domain)) &
+    call MOM_error(FATAL, " initialize_oda_incupd: Unable to open "//trim(filename))
+
+  call field_size(filename,h_var,siz,no_domain=.true.)
+  if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) &
+         call MOM_error(FATAL,"initialize_oda_incupd_file: Array size mismatch for oda data.")
+  nz_data = siz(3)
+  ! get h increments
+  allocate(hoda(isd:ied,jsd:jed,nz_data))
+  call MOM_read_data(filename, h_var   , hoda(:,:,:), G%Domain, scale=US%m_to_Z)
+  call initialize_oda_incupd( G, GV, US, param_file, oda_incupd_CSp, hoda, nz_data, restart_CS)
+  deallocate(hoda)
+
+  ! set-up T and S increments arrays
+  if (use_temperature) then
+    allocate(tmp_tr(isd:ied,jsd:jed,nz_data))
+    ! temperature inc. in array Inc(1)
+    call MOM_read_data(filename, tempinc_var, tmp_tr(:,:,:), G%Domain)
+    call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp)
+    ! salinity inc. in array Inc(2)
+    call MOM_read_data(filename, salinc_var, tmp_tr(:,:,:), G%Domain)
+    call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp)
+    deallocate(tmp_tr)
+  endif
+
+  ! set-up U and V increments arrays
+  if (uv_inc) then
+    filename = trim(inputdir)//trim(uv_inc_file)
+    call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename)
+    if (.not.file_exists(filename, G%Domain)) &
+            call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename))
+    allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data))
+    allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data))
+    tmp_u(:,:,:) = 0.0 ; tmp_v(:,:,:) = 0.0
+    call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T)
+    call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, oda_incupd_CSp)
+    deallocate(tmp_u,tmp_v)
+  endif
+
+  ! calculate increments if input are full fields
+  if (oda_inc) then ! input are increments
+    if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ")
+  else ! inputs are full fields
+    if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ")
+    call calc_oda_increments(h, tv, u, v, G, GV, US, oda_incupd_CSp)
+    if (save_inc) then
+      call output_oda_incupd_inc(Time, G, GV, param_file, oda_incupd_CSp, US)
+    endif
+  endif  ! not oda_inc
+
+end subroutine initialize_oda_incupd_file
+
+
 !> This subroutine sets the 4 bottom depths at velocity points to be the
 !! maximum of the adjacent depths.
 subroutine set_velocity_depth_max(G)