adjustment_initialization.F90

Go to the documentation of this file.

module adjustment_initialization
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type
use mom_io, only : close_file, fieldtype, file_exists
use mom_io, only : open_file, read_data, read_axis_data, single_file
use mom_io, only : write_field, slasher
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
use regrid_consts, only : coordinatemode, default_coordinate_mode
use regrid_consts, only : regridding_layer, regridding_zstar
use regrid_consts, only : regridding_rho, regridding_sigma

implicit none ; private

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

#include <MOM_memory.h>

! -----------------------------------------------------------------------------
! The following routines are visible to the outside world
! -----------------------------------------------------------------------------
public adjustment_initialize_thickness
public adjustment_initialize_temperature_salinity

! -----------------------------------------------------------------------------
! This module contains the following routines
! -----------------------------------------------------------------------------
contains

!------------------------------------------------------------------------------
!> Initialization of thicknesses.
!! This subroutine initializes the layer thicknesses to be uniform.
!------------------------------------------------------------------------------
subroutine adjustment_initialize_thickness ( h, G, GV, param_file, just_read_params)

  type(ocean_grid_type),   intent(in)  :: G           !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)  :: GV          !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized, in m.
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  real :: e0(szk_(g)+1)   ! The resting interface heights, in m, usually !
                          ! negative because it is positive upward.      !
  real :: eta1D(szk_(g)+1)! Interface height relative to the sea surface !
                          ! positive upward, in m.                       !
  real    :: x, y, yy, delta_S_strat, dSdz, delta_S, S_ref
  real    :: min_thickness, adjustment_width, adjustment_delta, adjustment_deltaS
  real    :: front_wave_amp, front_wave_length, front_wave_asym
  real    :: target_values(szk_(g)+1)
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=20) :: verticalCoordinate
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "adjustment_initialization"   ! This module's name.
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) &
    call mom_mesg("initialize_thickness_uniform: setting thickness")

  ! Parameters used by main model initialization
  if (.not.just_read) call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl,"S_REF",s_ref,fail_if_missing=.true.,do_not_log=.true.)
  call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness,'Minimum layer thickness', &
         units='m',default=1.0e-3, do_not_log=just_read)

  ! Parameters specific to this experiment configuration
  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
                 default=default_coordinate_mode, do_not_log=just_read)
  call get_param(param_file, mdl,"ADJUSTMENT_WIDTH",adjustment_width,     &
                 "Width of frontal zone",                                &
                 units="same as x,y", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"DELTA_S_STRAT",delta_s_strat,           &
                 "Top-to-bottom salinity difference of stratification",  &
                 units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"ADJUSTMENT_DELTAS",adjustment_deltas,   &
                 "Salinity difference across front",                     &
                 units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"FRONT_WAVE_AMP",front_wave_amp,         &
                 "Amplitude of trans-frontal wave perturbation",         &
                 units="same as x,y",default=0., do_not_log=just_read)
  call get_param(param_file, mdl,"FRONT_WAVE_LENGTH",front_wave_length,   &
                 "Wave-length of trans-frontal wave perturbation",       &
                 units="same as x,y",default=0., do_not_log=just_read)
  call get_param(param_file, mdl,"FRONT_WAVE_ASYM",front_wave_asym,       &
                 "Amplitude of frontal asymmetric perturbation",         &
                 default=0., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  ! WARNING: this routine specifies the interface heights so that the last layer
  !          is vanished, even at maximum depth. In order to have a uniform
  !          layer distribution, use this line of code within the loop:
  !          e0(k) = -G%max_depth * real(k-1) / real(nz)
  !          To obtain a thickness distribution where the last layer is
  !          vanished and the other thicknesses uniformly distributed, use:
  !          e0(k) = -G%max_depth * real(k-1) / real(nz-1)

  dsdz = -delta_s_strat/g%max_depth

  select case ( coordinatemode(verticalcoordinate) )

    case ( regridding_layer, regridding_rho )
      if (delta_s_strat.ne.0.) then
        adjustment_delta = adjustment_deltas / delta_s_strat * g%max_depth
        do k=1,nz+1
          e0(k) = adjustment_delta-(g%max_depth+2*adjustment_delta) * (real(k-1) / real(nz))
        enddo
      else
        adjustment_delta = 2.*g%max_depth
        do k=1,nz+1
          e0(k) = -(g%max_depth) * (real(k-1) / real(nz))
        enddo
      endif
      target_values(1)    = gv%Rlay(1)+0.5*(gv%Rlay(1)-gv%Rlay(2))
      target_values(nz+1) = gv%Rlay(nz)+0.5*(gv%Rlay(nz)-gv%Rlay(nz-1))
      do k = 2,nz
        target_values(k) = target_values(k-1) + ( gv%Rlay(nz) - gv%Rlay(1) ) / (nz-1)
      end do
      target_values = target_values - 1000.
      do j=js,je ; do i=is,ie
          if (front_wave_length.ne.0.) then
            y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
            yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / adjustment_width
            yy = min(1.0, yy); yy = max(-1.0, yy)
            yy = yy * 2. * acos( 0. )
            y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
          else
            y = 0.
          endif
          x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
          x = min(1.0, x); x = max(-1.0, x)
          x = x * acos( 0. )
          delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
          do k=2,nz
            if (dsdz.ne.0.) then
              eta1d(k) = ( target_values(k) - ( s_ref + delta_s ) ) / dsdz
            else
              eta1d(k) = e0(k) - (0.5*adjustment_delta) * sin( x )
            endif
            eta1d(k) = max( eta1d(k), -g%max_depth )
            eta1d(k) = min( eta1d(k), 0. )
          enddo
          eta1d(1)=0.; eta1d(nz+1)=-g%max_depth
          do k=nz,1,-1
            if (eta1d(k) > 0.) then
              eta1d(k) = max( eta1d(k+1) + min_thickness, 0. )
              h(i,j,k) = max( eta1d(k) - eta1d(k+1), min_thickness )
            elseif (eta1d(k) <= (eta1d(k+1) + min_thickness)) then
              eta1d(k) = eta1d(k+1) + min_thickness
              h(i,j,k) = min_thickness
            else
              h(i,j,k) = eta1d(k) - eta1d(k+1)
            endif
          enddo
      enddo ; enddo

    case ( regridding_zstar, regridding_sigma )
      do k=1,nz+1
        eta1d(k) = -(g%max_depth) * (real(k-1) / real(nz))
        eta1d(k) = max(min(eta1d(k),0.),-g%max_depth)
      enddo
      do j=js,je ; do i=is,ie
        do k=nz,1,-1
            h(i,j,k) = eta1d(k) - eta1d(k+1)
        enddo
      enddo ; enddo

    case default
      call mom_error(fatal,"adjustment_initialize_thickness: "// &
      "Unrecognized i.c. setup - set ADJUSTMENT_IC")

  end select

end subroutine adjustment_initialize_thickness


!------------------------------------------------------------------------------
!> Initialization of temperature and salinity.
!------------------------------------------------------------------------------
subroutine adjustment_initialize_temperature_salinity ( T, S, h, G, param_file, &
                                                    eqn_of_state, just_read_params)
  type(ocean_grid_type),   intent(in) :: G                    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< The temperature that is being initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< The salinity that is being initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< The model thickness.
  type(param_file_type),   intent(in) :: param_file           !< A structure indicating the
                                         !! open file to parse for model parameter values.
  type(eos_type),                 pointer     :: eqn_of_state !< Equation of state.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  integer   :: i, j, k, is, ie, js, je, nz
  real      :: x, y, yy
  integer   :: index_bay_z
  real      :: S_ref, T_ref         ! Reference salinity and temerature within
                                    ! surface layer
  real      :: S_range, T_range     ! Range of salinities and temperatures over the
                                    ! vertical
  real      :: xi0, xi1, dSdz, delta_S, delta_S_strat
  real      :: adjustment_width, adjustment_deltaS
  real       :: front_wave_amp, front_wave_length, front_wave_asym
  real      :: eta1d(szk_(g)+1)
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=20) :: verticalCoordinate

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  ! Parameters used by main model initialization
  call get_param(param_file, mdl,"S_REF",s_ref,'Reference salinity', units='1e-3', &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"T_REF",t_ref,'Reference temperature', units='C', &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"S_RANGE",s_range,'Initial salinity range', units='1e-3', &
                 default=2.0, do_not_log=just_read)
  call get_param(param_file, mdl,"T_RANGE",t_range,'Initial temperature range', units='C', &
                 default=0.0, do_not_log=just_read)
  ! Parameters specific to this experiment configuration BUT logged in previous s/r
  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
                 default=default_coordinate_mode, do_not_log=just_read)
  call get_param(param_file, mdl,"ADJUSTMENT_WIDTH", adjustment_width, &
                 fail_if_missing=.not.just_read, do_not_log=.true.)
  call get_param(param_file, mdl,"ADJUSTMENT_DELTAS", adjustment_deltas, &
                 fail_if_missing=.not.just_read, do_not_log=.true.)
  call get_param(param_file, mdl,"DELTA_S_STRAT", delta_s_strat, &
                 fail_if_missing=.not.just_read, do_not_log=.true.)
  call get_param(param_file, mdl,"FRONT_WAVE_AMP", front_wave_amp, default=0., &
                 do_not_log=.true.)
  call get_param(param_file, mdl,"FRONT_WAVE_LENGTH",front_wave_length, &
                 default=0.,do_not_log=.true.)
  call get_param(param_file, mdl,"FRONT_WAVE_ASYM", front_wave_asym, default=0., &
                 do_not_log=.true.)

  if (just_read) return ! All run-time parameters have been read, so return.

  t(:,:,:) = 0.0
  s(:,:,:) = 0.0

  ! Linear salinity profile
  select case ( coordinatemode(verticalcoordinate) )

    case ( regridding_zstar, regridding_sigma )
      dsdz = -delta_s_strat/g%max_depth
      do j=js,je ; do i=is,ie
          eta1d(nz+1)=-g%bathyT(i,j)
          do k=nz,1,-1
            eta1d(k)=eta1d(k+1)+h(i,j,k)
          enddo
          if (front_wave_length.ne.0.) then
            y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
            yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / front_wave_length
            yy = min(1.0, yy); yy = max(-1.0, yy)
            yy = yy * 2. * acos( 0. )
            y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
          else
            y = 0.
          endif
          x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
          x = min(1.0, x); x = max(-1.0, x)
          x = x * acos( 0. )
          delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
          do k=1,nz
            s(i,j,k) = s_ref + delta_s + 0.5 * ( eta1d(k)+eta1d(k+1) ) * dsdz
            x = abs(s(i,j,k) - 0.5*real(nz-1)/real(nz)*S_range)/s_range*real(2*nz)
            x = 1.-min(1., x)
            t(i,j,k) = x
         enddo
   !     x=sum(T(i,j,:)*h(i,j,:))
   !     T(i,j,:)=T(i,j,:)/x*(G%max_depth*1.5/real(nz))
      enddo ; enddo

    case ( regridding_layer, regridding_rho )
      do k = 1,nz
        s(:,:,k) = s_ref + s_range * ( (real(k)-0.5) / real( nz ) )
   !    x = abs(S(1,1,k) - 0.5*real(nz-1)/real(nz)*S_range)/S_range*real(2*nz)
   !    x = 1.-min(1., x)
   !    T(:,:,k) = x
      end do

    case default
      call mom_error(fatal,"adjustment_initialize_temperature_salinity: "// &
      "Unrecognized i.c. setup - set ADJUSTMENT_IC")

  end select

end subroutine adjustment_initialize_temperature_salinity

!> \namespace adjustment_initialization
!!
!! The module configures the model for the geostrophic adjustment
!! test case.
end module adjustment_initialization