phillips_initialization Module Reference

Detailed Description

Initialization for the "Phillips" channel configuration.

By Robert Hallberg, April 1994 - June 2002

This subroutine initializes the fields for the simulations. The one argument passed to initialize, Time, is set to the current time of the simulation. The fields which are initialized here are: u - Zonal velocity [m s-1]. v - Meridional velocity [m s-1]. h - Layer thickness in m. (Must be positive.) D - Basin depth in m. (Must be positive.) f - The Coriolis parameter [s-1]. g - The reduced gravity at each interface [m s-2] Rlay - Layer potential density (coordinate variable) [kg m-3]. If ENABLE_THERMODYNAMICS is defined: T - Temperature [degC]. S - Salinity [ppt]. If SPONGE is defined: A series of subroutine calls are made to set up the damping rates and reference profiles for all variables that are damped in the sponge. Any user provided tracer code is also first linked through this subroutine.

Forcing-related fields (taux, tauy, buoy, ustar, etc.) are set in MOM_surface_forcing.F90.

These variables are all set in the set of subroutines (in this file) Phillips_initialize_thickness, Phillips_initialize_velocity, Phillips_initialize_topography and Phillips_initialize_sponges that seet up fields that are specific to the Phillips instability test case.

Functions/Subroutines
subroutine, public	phillips_initialize_thickness (h, G, GV, US, param_file, just_read_params)
	Initialize the thickness field for the Phillips model test case. More...
subroutine, public	phillips_initialize_velocity (u, v, G, GV, US, param_file, just_read_params)
	Initialize the velocity fields for the Phillips model test case. More...
subroutine, public	phillips_initialize_sponges (G, GV, US, tv, param_file, CSp, h)
	Sets up the the inverse restoration time (Idamp), and the values towards which the interface heights and an arbitrary number of tracers should be restored within each sponge for the Phillips model test case. More...
real function	sech (x)
	sech calculates the hyperbolic secant. More...
subroutine, public	phillips_initialize_topography (D, G, param_file, max_depth, US)
	Initialize topography. More...

Function/Subroutine Documentation

phillips_initialize_sponges()

subroutine, public phillips_initialization::phillips_initialize_sponges	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(thermo_var_ptrs), intent(in)	tv,
		type(param_file_type), intent(in)	param_file,
		type(sponge_cs), pointer	CSp,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	h
	)

Sets up the the inverse restoration time (Idamp), and the values towards which the interface heights and an arbitrary number of tracers should be restored within each sponge for the Phillips model test case.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	tv	A structure containing pointers to any available thermodynamic fields, potential temperature and salinity or mixed layer density. Absent fields have NULL ptrs.
[in]	param_file	A structure indicating the open file to parse for model parameter values.
	csp	A pointer that is set to point to the control structure for the sponge module.
[in]	h	Thickness field [H ~> m or kg m-2].

Definition at line 198 of file Phillips_initialization.F90.

  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure.
  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
  type(unit_scale_type), intent(in) :: US   !< A dimensional unit scaling type
  type(thermo_var_ptrs), intent(in) :: tv   !< A structure containing pointers
                                            !! to any available thermodynamic
                                            !! fields, potential temperature and
                                            !! salinity or mixed layer density.
                                            !! Absent fields have NULL ptrs.
  type(param_file_type), intent(in) :: param_file !< A structure indicating the
                                            !! open file to parse for model
                                            !! parameter values.
  type(sponge_CS),   pointer    :: CSp      !< A pointer that is set to point to
                                            !! the control structure for the
                                            !! sponge module.
  real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Thickness field [H ~> m or kg m-2].
 
  ! Local variables
  real :: eta0(SZK_(G)+1)   ! The 1-d nominal positions of the interfaces.
  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta [Z ~> m].
  real :: temp(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for other variables.
  real :: Idamp(SZI_(G),SZJ_(G))    ! The inverse damping rate [s-1].
  real :: eta_im(SZJ_(G),SZK_(G)+1) ! A temporary array for zonal-mean eta [Z ~> m].
  real :: Idamp_im(SZJ_(G))         ! The inverse zonal-mean damping rate [s-1].
  real :: damp_rate    ! The inverse zonal-mean damping rate [s-1].
  real :: jet_width    ! The width of the zonal mean jet, in km.
  real :: jet_height   ! The interface height scale associated with the zonal-mean jet [Z ~> m].
  real :: y_2          ! The y-position relative to the channel center, in km.
  real :: half_strat   ! The fractional depth where the straficiation is centered [nondim].
  real :: half_depth   ! The depth where the stratification is centered [Z ~> m].
  character(len=40)  :: mdl = "Phillips_initialize_sponges" ! This subroutine's name.
 
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
  logical, save :: first_call = .true.
 
  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
 
  eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
  eta_im(:,:) = 0.0 ; idamp_im(:) = 0.0
 
  if (first_call) call log_version(param_file, mdl, version)
  first_call = .false.
  call get_param(param_file, mdl, "HALF_STRAT_DEPTH", half_strat, &
                 "The fractional depth where the stratificaiton is centered.", &
                 units="nondim", default = 0.5)
  call get_param(param_file, mdl, "SPONGE_RATE", damp_rate, &
                 "The rate at which the zonal-mean sponges damp.", units="s-1", &
                 default = 1.0/(10.0*86400.0))
 
  call get_param(param_file, mdl, "JET_WIDTH", jet_width, &
                 "The width of the zonal-mean jet.", units="km", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
                 "The interface height scale associated with the "//&
                 "zonal-mean jet.", units="m", scale=us%m_to_Z, &
                 fail_if_missing=.true.)
 
  half_depth = g%max_depth*half_strat
  eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
  do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ; enddo
  do k=2+nz/2,nz+1
    eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
  enddo
 
  do j=js,je
    idamp_im(j) = damp_rate
    eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
  enddo
  do k=2,nz ; do j=js,je
    y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
    eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
!         jet_height * atan(y_2 / jet_width)
    if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
    if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
  enddo ; enddo
 
  call initialize_sponge(idamp, eta, g, param_file, csp, gv, idamp_im, eta_im)
 

phillips_initialize_thickness()

subroutine, public phillips_initialization::phillips_initialize_thickness	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out)	h,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		logical, intent(in), optional	just_read_params
	)

Initialize the thickness field for the Phillips model test case.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[out]	h	The thickness that is being initialized [H ~> m or kg m-2]
[in]	param_file	A structure indicating the open file to parse for model parameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 39 of file Phillips_initialization.F90.

  type(ocean_grid_type),   intent(in)  :: 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
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h          !< The thickness that is being initialized [H ~> m or kg m-2]
  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 :: eta0(SZK_(G)+1)   ! The 1-d nominal positions of the interfaces [Z ~> m]
  real :: eta_im(SZJ_(G),SZK_(G)+1) ! A temporary array for zonal-mean eta [Z ~> m]
  real :: eta1D(SZK_(G)+1)  ! Interface height relative to the sea surface, positive upward [Z ~> m]
  real :: jet_width         ! The width of the zonal-mean jet [km]
  real :: jet_height        ! The interface height scale associated with the zonal-mean jet [Z ~> m]
  real :: y_2
  real :: half_strat      ! The fractional depth where the stratification is centered [nondim]
  real :: half_depth      ! The depth where the stratification is centered [Z ~> m]
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "Phillips_initialize_thickness" ! This subroutine's name.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
 
  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
 
  eta_im(:,:) = 0.0
 
  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
 
  if (.not.just_read) call log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "HALF_STRAT_DEPTH", half_strat, &
                 "The fractional depth where the stratification is centered.", &
                 units="nondim", default = 0.5, do_not_log=just_read)
  call get_param(param_file, mdl, "JET_WIDTH", jet_width, &
                 "The width of the zonal-mean jet.", units="km", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
                 "The interface height scale associated with the "//&
                 "zonal-mean jet.", units="m", scale=us%m_to_Z, &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
 
  if (just_read) return ! All run-time parameters have been read, so return.
 
  half_depth = g%max_depth*half_strat
  eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
  do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ; enddo
  do k=2+nz/2,nz+1
    eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
  enddo
 
  do j=js,je
    eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
  enddo
  do k=2,nz ; do j=js,je
    y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
    eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
                ! or  ... + jet_height * atan(y_2 / jet_width)
    if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
    if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
  enddo ; enddo
 
  do j=js,je ; do i=is,ie
    !   This sets the initial thickness in [H ~> m or kg m-2] of the layers.  The
    ! thicknesses are set to insure that: 1. each layer is at least an Angstrom thick, and
    ! 2. the interfaces are where they should be based on the resting depths and interface
    !    height perturbations, as long at this doesn't interfere with 1.
    eta1d(nz+1) = -g%bathyT(i,j)
    do k=nz,1,-1
      eta1d(k) = eta_im(j,k)
      if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
        eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
        h(i,j,k) = gv%Angstrom_H
      else
        h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
      endif
    enddo
  enddo ; enddo
 

phillips_initialize_topography()

subroutine, public phillips_initialization::phillips_initialize_topography	(	real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)	D,
		type(dyn_horgrid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		real, intent(in)	max_depth,
		type(unit_scale_type), intent(in), optional	US
	)

Initialize topography.

Parameters

[in]	g	The dynamic horizontal grid type
[out]	d	Ocean bottom depth in m or Z if US is present
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum model depth in the units of D
[in]	us	A dimensional unit scaling type

Definition at line 293 of file Phillips_initialization.F90.

  type(dyn_horgrid_type),          intent(in)  :: G !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                   intent(out) :: D !< Ocean bottom depth in m or Z if US is present
  type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
  real,                            intent(in)  :: max_depth !< Maximum model depth in the units of D
  type(unit_scale_type), optional, intent(in)  :: US !< A dimensional unit scaling type
 
  ! Local variables
  real :: m_to_Z  ! A dimensional rescaling factor.
  real :: PI, Htop, Wtop, Ltop, offset, dist
  real :: x1, x2, x3, x4, y1, y2
  integer :: i,j,is,ie,js,je
  character(len=40)  :: mdl = "Phillips_initialize_topography" ! This subroutine's name.
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  pi = 4.0*atan(1.0)
  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
 
  call get_param(param_file, mdl, "PHILLIPS_HTOP", htop, &
                 "The maximum height of the topography.", units="m", scale=m_to_z, &
                 fail_if_missing=.true.)
! Htop=0.375*max_depth     ! max height of topog. above max_depth
  wtop = 0.5*g%len_lat     ! meridional width of drake and mount
  ltop = 0.25*g%len_lon    ! zonal width of topographic features
  offset = 0.1*g%len_lat   ! meridional offset from center
  dist = 0.333*g%len_lon   ! distance between drake and mount
                           ! should be longer than Ltop/2
 
  y1=g%south_lat+0.5*g%len_lat+offset-0.5*wtop; y2=y1+wtop
  x1=g%west_lon+0.1*g%len_lon; x2=x1+ltop; x3=x1+dist; x4=x3+3.0/2.0*ltop
 
  do i=is,ie ; do j=js,je
    d(i,j)=0.0
    if (g%geoLonT(i,j)>x1 .and. g%geoLonT(i,j)<x2) then
      d(i,j) = htop*sin(pi*(g%geoLonT(i,j)-x1)/(x2-x1))**2
      if (g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2) then
         d(i,j) = d(i,j)*(1-sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2)
      endif
    elseif (g%geoLonT(i,j)>x3 .and. g%geoLonT(i,j)<x4 .and. &
             g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2) then
      d(i,j) = 2.0/3.0*htop*sin(pi*(g%geoLonT(i,j)-x3)/(x4-x3))**2 &
                   *sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2
    endif
    d(i,j) = max_depth - d(i,j)
  enddo ; enddo
 

phillips_initialize_velocity()

subroutine, public phillips_initialization::phillips_initialize_velocity	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)	v,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		logical, intent(in), optional	just_read_params
	)

Initialize the velocity fields for the Phillips model test case.

Parameters

[in]	g	Grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[out]	u	i-component of velocity [m s-1]
[out]	v	j-component of velocity [m s-1]
[in]	param_file	A structure indicating the open file to parse for modelparameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 121 of file Phillips_initialization.F90.

  type(ocean_grid_type),   intent(in)  :: G  !< Grid structure
  type(verticalGrid_type), intent(in)  :: GV !< Vertical grid structure
  type(unit_scale_type),   intent(in)  :: US !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out) :: u  !< i-component of velocity [m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out) :: v  !< j-component of velocity [m s-1]
  type(param_file_type),   intent(in)  :: param_file !< A structure indicating the open file to
                                                     !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                     !! only read parameters without changing h.
 
  real :: jet_width, jet_height, x_2, y_2
  real :: velocity_amplitude, pi
  integer :: i, j, k, is, ie, js, je, nz, m
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "Phillips_initialize_velocity" ! This subroutine's name.
  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 log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "VELOCITY_IC_PERTURB_AMP", velocity_amplitude, &
                 "The magnitude of the initial velocity perturbation.", &
                 units="m s-1", default=0.001, do_not_log=just_read)
  call get_param(param_file, mdl, "JET_WIDTH", jet_width, &
                 "The width of the zonal-mean jet.", units="km", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
                 "The interface height scale associated with the "//&
                 "zonal-mean jet.", units="m", scale=us%m_to_Z, &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
 
  if (just_read) return ! All run-time parameters have been read, so return.
 
  u(:,:,:) = 0.0
  v(:,:,:) = 0.0
 
  pi = 4.0*atan(1.0)
 
  ! Use thermal wind shear to give a geostrophically balanced flow.
  do k=nz-1,1 ; do j=js,je ; do i=is-1,ie
    y_2 = g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat
! This uses d/d y_2 atan(y_2 / jet_width)
!    u(I,j,k) = u(I,j,k+1) + (1e-3 * jet_height / &
!           (jet_width * (1.0 + (y_2 / jet_width)**2))) * &
!           (2.0 * US%L_to_m**2*US%s_to_T**2*GV%g_prime(K+1) * US%T_to_s / (G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)))
! This uses d/d y_2 tanh(y_2 / jet_width)
    u(i,j,k) = u(i,j,k+1) + (1e-3 * (jet_height / jet_width) * &
           (sech(y_2 / jet_width))**2 ) * &
           (2.0 * us%L_to_m**2*us%s_to_T**2*gv%g_prime(k+1) * us%T_to_s / (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)))
  enddo ; enddo ; enddo
 
  do k=1,nz ; do j=js,je ; do i=is-1,ie
    y_2 = (g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat) / g%len_lat
    x_2 = (g%geoLonCu(i,j) - g%west_lon - 0.5*g%len_lon) / g%len_lon
    if (g%geoLonCu(i,j) == g%west_lon) then
      ! This modification is required so that the perturbations are identical for
      ! symmetric and non-symmetric memory.  It is exactly equivalent to
      ! taking the longitude at the eastern edge of the domain, so that x_2 ~= 0.5.
      x_2 = ((g%west_lon + g%len_lon*real(g%ieg-(g%isg-1))/real(g%Domain%niglobal)) - &
             g%west_lon - 0.5*g%len_lon) / g%len_lon
    endif
    u(i,j,k) = u(i,j,k) + velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
           (0.5 - abs(2.0*x_2) + 0.1*abs(cos(10.0*pi*x_2)) - abs(sin(5.0*pi*y_2)))
    do m=1,10
      u(i,j,k) = u(i,j,k) + 0.2*velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
            cos(2.0*m*pi*x_2 + 2*m) * cos(6.0*pi*y_2)
    enddo
  enddo ; enddo ; enddo
 

sech()

real function phillips_initialization::sech ( real, intent(in)  x )

private

sech calculates the hyperbolic secant.

Parameters

[in]

x

Input value.

Returns

Result.

Definition at line 280 of file Phillips_initialization.F90.

  real, intent(in) :: x    !< Input value.
  real             :: sech !< Result.
 
  ! This is here to prevent overflows or underflows.
  if (abs(x) > 228.) then
    sech = 0.0
  else
    sech = 2.0 / (exp(x) + exp(-x))
  endif