MOM_tracer_initialization_from_Z.F90

!> Initializes hydrography from z-coordinate climatology files
module mom_tracer_initialization_from_z

! This file is part of MOM6. See LICENSE.md for the license.

use mom_debugging, only : hchksum
use mom_coms, only : max_across_pes, min_across_pes
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only :  clock_routine, clock_loop
use mom_density_integrals, only : int_specific_vol_dp
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 : 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
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type, ispointincell
use mom_horizontal_regridding, only : mystats, horiz_interp_and_extrap_tracer
use mom_regridding, only : regridding_cs
use mom_remapping, only : remapping_cs, initialize_remapping
use mom_remapping, only : remapping_core_h
use mom_string_functions, only : uppercase
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type, setverticalgridaxes
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
use mom_ale, only : ale_remap_scalar

implicit none ; private

#include <MOM_memory.h>

public :: mom_initialize_tracer_from_z

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

character(len=40)  :: mdl = "MOM_tracer_initialization_from_Z" !< This module's name.

contains

!> Initializes a tracer from a z-space data file.
subroutine mom_initialize_tracer_from_z(h, tr, G, GV, US, PF, src_file, src_var_nam, &
                          src_var_unit_conversion, src_var_record, homogenize, &
                          useALEremapping, remappingScheme, src_var_gridspec )
  type(ocean_grid_type),      intent(inout) :: G   !< Ocean grid structure.
  type(verticalgrid_type),    intent(in)    :: GV  !< Ocean vertical grid structure.
  type(unit_scale_type),      intent(in)    :: US  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                              intent(in)    :: h   !< Layer thickness [H ~> m or kg m-2].
  real, dimension(:,:,:),     pointer       :: tr  !< Pointer to array to be initialized
  type(param_file_type),      intent(in)    :: PF  !< parameter file
  character(len=*),           intent(in)    :: src_file !< source filename
  character(len=*),           intent(in)    :: src_var_nam !< variable name in file
  real,             optional, intent(in)    :: src_var_unit_conversion !< optional multiplicative unit conversion
  integer,          optional, intent(in)    :: src_var_record  !< record to read for multiple time-level files
  logical,          optional, intent(in)    :: homogenize !< optionally homogenize to mean value
  logical,          optional, intent(in)    :: useALEremapping !< to remap or not (optional)
  character(len=*), optional, intent(in)    :: remappingScheme !< remapping scheme to use.
  character(len=*), optional, intent(in)    :: src_var_gridspec !< Source variable name in a gridspec file.
                                                                !! This is not implemented yet.
  ! Local variables
  real :: land_fill = 0.0
  character(len=200) :: inputdir ! The directory where NetCDF input files are.
  character(len=200) :: mesg
  real               :: convert
  integer            :: recnum
  character(len=10)  :: remapScheme
  logical            :: homog,useALE

  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_initialize_tracers_from_Z" ! This module's name.

  integer :: is, ie, js, je, nz ! compute domain indices
  integer :: isd, ied, jsd, jed ! data domain indices
  integer :: i, j, k, kd
  real, allocatable, dimension(:,:,:), target :: tr_z, mask_z
  real, allocatable, dimension(:), target :: z_edges_in, z_in

  ! Local variables for ALE remapping
  real, dimension(:,:,:), allocatable :: hSrc ! Source thicknesses [H ~> m or kg m-2].
  real, dimension(:), allocatable :: h1 ! A 1-d column of source thicknesses [Z ~> m].
  real :: zTopOfCell, zBottomOfCell, z_bathy  ! Heights [Z ~> m].
  type(remapping_cs) :: remapCS ! Remapping parameters and work arrays

  real :: missing_value
  integer :: nPoints
  integer :: id_clock_routine, id_clock_ALE
  logical :: answers_2018, default_2018_answers, hor_regrid_answers_2018
  logical :: reentrant_x, tripolar_n

  id_clock_routine = cpu_clock_id('(Initialize tracer from Z)', grain=clock_routine)
  id_clock_ale = cpu_clock_id('(Initialize tracer from Z) ALE', grain=clock_loop)

  call cpu_clock_begin(id_clock_routine)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homog, &
                 "If True, then horizontally homogenize the interpolated "//&
                 "initial conditions.", default=.false.)
  call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", useale, &
                 "If True, then remap straight to model coordinate from file.",&
                 default=.true.)
  call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remapscheme, &
                 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING is True.", &
                 default="PLM")
  call get_param(pf, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  if (useale) then
    call get_param(pf, mdl, "REMAPPING_2018_ANSWERS", answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
  endif
  call get_param(pf, mdl, "HOR_REGRID_2018_ANSWERS", hor_regrid_answers_2018, &
                 "If true, use the order of arithmetic for horizonal regridding that recovers "//&
                 "the answers from the end of 2018.  Otherwise, use rotationally symmetric "//&
                 "forms of the same expressions.", default=default_2018_answers)

  ! These are model grid properties, but being applied to the data grid for now.
  ! need to revisit this (mjh)
  reentrant_x = .false. ;  call get_param(pf, mdl, "REENTRANT_X", reentrant_x,default=.true.)
  tripolar_n = .false. ;  call get_param(pf, mdl, "TRIPOLAR_N", tripolar_n, default=.false.)

  if (PRESENT(homogenize)) homog=homogenize
  if (PRESENT(usealeremapping)) useale=usealeremapping
  if (PRESENT(remappingscheme)) remapscheme=remappingscheme
  recnum=1
  if (PRESENT(src_var_record)) recnum = src_var_record
  convert=1.0
  if (PRESENT(src_var_unit_conversion)) convert = src_var_unit_conversion

  call horiz_interp_and_extrap_tracer(src_file, src_var_nam, convert, recnum, &
       g, tr_z, mask_z, z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, &
       homog, m_to_z=us%m_to_Z, answers_2018=hor_regrid_answers_2018)

  kd = size(z_edges_in,1)-1
  call pass_var(tr_z,g%Domain)
  call pass_var(mask_z,g%Domain)

! Done with horizontal interpolation.
! Now remap to model coordinates
  if (useale) then
    call cpu_clock_begin(id_clock_ale)
    ! First we reserve a work space for reconstructions of the source data
    allocate( h1(kd) )
    allocate( hsrc(isd:ied,jsd:jed,kd) )
    ! Set parameters for reconstructions
    call initialize_remapping( remapcs, remapscheme, boundary_extrapolation=.false., answers_2018=answers_2018 )
    ! Next we initialize the regridding package so that it knows about the target grid

    do j = js, je ; do i = is, ie
      if (g%mask2dT(i,j)>0.) then
        ! Build the source grid
        ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0
        z_bathy = g%bathyT(i,j)
        do k = 1, kd
          if (mask_z(i,j,k) > 0.) then
            zbottomofcell = -min( z_edges_in(k+1), z_bathy )
          elseif (k>1) then
            zbottomofcell = -z_bathy
          endif
          h1(k) = ztopofcell - zbottomofcell
          if (h1(k)>0.) npoints = npoints + 1
          ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
        enddo
        h1(kd) = h1(kd) + ( ztopofcell + z_bathy ) ! In case data is deeper than model
      else
        tr(i,j,:) = 0.
      endif ! mask2dT
      hsrc(i,j,:) = gv%Z_to_H * h1(:)
    enddo ; enddo

    call ale_remap_scalar(remapcs, g, gv, kd, hsrc, tr_z, h, tr, all_cells=.false., answers_2018=answers_2018 )

    deallocate( hsrc )
    deallocate( h1 )

    do k=1,nz
      call mystats(tr(:,:,k), missing_value, is, ie, js, je, k, 'Tracer from ALE()')
    enddo
    call cpu_clock_end(id_clock_ale)
  endif ! useALEremapping

! Fill land values
  do k=1,nz ; do j=js,je ; do i=is,ie
    if (tr(i,j,k) == missing_value) then
      tr(i,j,k)=land_fill
    endif
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')
  call cpu_clock_end(id_clock_routine)

end subroutine mom_initialize_tracer_from_z

end module mom_tracer_initialization_from_z