phillips_initialization Module Reference

Detailed Description

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 in m s-1. * v - Meridional velocity in m s-1. * h - Layer thickness in m. (Must be positive.) * D - Basin depth in m. (Must be positive.) * f - The Coriolis parameter, in s-1. * g - The reduced gravity at each interface, in m s-2. * Rlay - Layer potential density (coordinate variable) in kg m-3. * If ENABLE_THERMODYNAMICS is defined: * T - Temperature in C. * S - Salinity in psu. * 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. *
• Macros written all in capital letters are defined in MOM_memory.h. *
• A small fragment of the grid is shown below: *
• j+1 x ^ x ^ x At x: q, f * j+1 > o > o > At ^: v, tauy * j x ^ x ^ x At >: u, taux * j > o > o > At o: h, D, buoy, tr, T, S, ustar * j-1 x ^ x ^ x * i-1 i i+1 At x & ^: * i i+1 At > & o: *
• The boundaries always run through q grid points (x). *.

Functions/Subroutines
subroutine, public	phillips_initialize_thickness (h, G, GV, param_file, just_read_params)
	Initialize thickness field. More...
subroutine, public	phillips_initialize_velocity (u, v, G, GV, param_file, just_read_params)
	Initialize velocity fields. More...
subroutine, public	phillips_initialize_sponges (G, use_temperature, tv, param_file, CSp, h)
	Sets up the the inverse restoration time (Idamp), and. More...
real function	sech (x)
	sech calculates the hyperbolic secant. More...
subroutine, public	phillips_initialize_topography (D, G, param_file, max_depth)
	Initialize topography. More...

Function/Subroutine Documentation

phillips_initialize_sponges()

subroutine, public phillips_initialization::phillips_initialize_sponges	(	type(ocean_grid_type), intent(in)	G,
		logical, intent(in)	use_temperature,
		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, g %ke), intent(in)	h
	)

Sets up the the inverse restoration time (Idamp), and.

Parameters

[in]	g	The ocean's grid structure.
[in]	use_temperature	Switch for temperature.
[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.

Definition at line 209 of file Phillips_initialization.F90.

References mom_sponge::initialize_sponge().

Referenced by mom_state_initialization::mom_initialize_state().

  type(ocean_grid_type), intent(in) :: g    !< The ocean's grid structure.
  logical, intent(in) :: use_temperature    !< Switch for temperature.
  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_(G)) :: h !< Thickness field.

  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, 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, in s-1.
  real :: eta_im(szj_(g),szk_(g)+1) ! A temporary array for zonal-mean eta, m.
  real :: idamp_im(szj_(g))         ! The inverse zonal-mean damping rate, in s-1.
  real :: damp_rate, jet_width, jet_height, y_2
  real :: half_strat, half_depth
  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 maximum depth of the ocean.", 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 \n"//&
                 "zonal-mean jet.", units="m", &
                 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, idamp_im, eta_im)

Here is the call graph for this function:

Here is the caller graph for this function:

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(param_file_type), intent(in)	param_file,
		logical, intent(in), optional	just_read_params
	)

Initialize thickness field.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[out]	h	The thickness that is being initialized, in m.
[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 52 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.
  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 :: eta0(szk_(g)+1)   ! The 1-d nominal positions of the interfaces.
  real :: eta_im(szj_(g),szk_(g)+1) ! A temporary array for zonal-mean eta, m.
  real :: eta1d(szk_(g)+1)  ! Interface height relative to the sea surface
                            ! positive upward, in m.
  real :: damp_rate, jet_width, jet_height, y_2
  real :: half_strat, half_depth
  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 maximum depth of the ocean.", 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 \n"//&
                 "zonal-mean jet.", units="m", &
                 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)
!         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 m) 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) = -1.0*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_z
      else
        h(i,j,k) = 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
	)

Initialize topography.

Parameters

[in]	g	The dynamic horizontal grid type
[out]	d	Ocean bottom depth in m
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum depth of model in m

Definition at line 299 of file Phillips_initialization.F90.

Referenced by mom_fixed_initialization::mom_initialize_topography().

  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
  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m

  real :: pi, htop, wtop, ltop, offset, dist, &
          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)

  call get_param(param_file, mdl, "PHILLIPS_HTOP", htop,             &
                 "The maximum height of the topography.", units="m", &
                 fail_if_missing=.true.)
! Htop=0.375*G%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)
       end if
     else if (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
     end if
     d(i,j)=max_depth-d(i,j)
  enddo; enddo

Here is the caller graph for this function:

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(param_file_type), intent(in)	param_file,
		logical, intent(in), optional	just_read_params
	)

Initialize velocity fields.

Parameters

[in]	g	Grid structure
[in]	gv	Vertical grid structure
[out]	u	i-component of velocity [m/s]
[out]	v	j-component of velocity [m/s]
[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 133 of file Phillips_initialization.F90.

References sech().

  type(ocean_grid_type),   intent(in)  :: g  !< Grid structure
  type(verticalgrid_type), intent(in)  :: gv !< Vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out) :: u  !< i-component of velocity [m/s]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out) :: v  !< j-component of velocity [m/s]
  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 :: damp_rate, 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 \n"//&
                 "zonal-mean jet.", units="m", &
                 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 * GV%g_prime(K+1) / (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 * gv%g_prime(k+1) / (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

Here is the call graph for this function:

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 286 of file Phillips_initialization.F90.

Referenced by phillips_initialize_velocity().

  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

Here is the caller graph for this function: