scm_idealized_hurricane Module Reference

Detailed Description

Initial conditions and forcing for the single column model (SCM) idealized hurricane example.

Data Types
type	scm_idealized_hurricane_cs
	Container for parameters describing idealized wind structure. More...

Functions/Subroutines
subroutine, public	scm_idealized_hurricane_ts_init (T, S, h, G, GV, param_file, just_read_params)
	Initializes temperature and salinity for the SCM idealized hurricane example. More...
subroutine, public	scm_idealized_hurricane_wind_init (Time, G, param_file, CS)
	Initializes wind profile for the SCM idealized hurricane example. More...
subroutine, public	scm_idealized_hurricane_wind_forcing (state, fluxes, day, G, CS)

Variables
character(len=40)	mdl = "SCM_idealized_hurricane"

Function/Subroutine Documentation

scm_idealized_hurricane_ts_init()

subroutine, public scm_idealized_hurricane::scm_idealized_hurricane_ts_init	(	real, dimension(szi_(g),szj_(g), szk_(g)), intent(out)	T,
		real, dimension(szi_(g),szj_(g), szk_(g)), intent(out)	S,
		real, dimension(szi_(g),szj_(g), szk_(g)), intent(in)	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
	)

Initializes temperature and salinity for the SCM idealized hurricane example.

Parameters

[in]	g	Grid structure
[in]	gv	Vertical grid structure
[out]	t	Potential temperature (degC)
[out]	s	Salinity (psu)
[in]	h	Layer thickness (m or Pa)
[in]	param_file	Input parameter structure
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 48 of file SCM_idealized_hurricane.F90.

References mdl.

Referenced by mom_state_initialization::mom_initialize_state().

  type(ocean_grid_type),                     intent(in)  :: g !< Grid structure
  type(verticalgrid_type),                   intent(in)  :: gv !< Vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< Potential temperature (degC)
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< Salinity (psu)
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness (m or Pa)
  type(param_file_type),                     intent(in)  :: param_file !< Input parameter structure
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
  ! Local variables
  real :: eta(szk_(g)+1) ! The 1-d nominal positions of the interfaces.
  real :: s_ref, sst_ref, dtdz, mld
  real :: zc
  logical :: just_read    ! If true, just read parameters but set nothing.
  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

  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,"SCM_S_REF",s_ref, &
                 'Reference salinity', units='1e-3',default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_SST_REF",sst_ref, &
                 'Reference surface temperature', units='C', &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_DTDZ",dtdz,                         &
                 'Initial temperature stratification below mixed layer', &
                 units='C/m', fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_MLD",mld, &
                 'Initial mixed layer depth', 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.

  do j=js,je ; do i=is,ie
    eta(1) = 0. ! Reference to surface
    do k=1,nz
      eta(k+1) = eta(k) - h(i,j,k)*gv%H_to_m ! Interface below layer (in m)
      zc = 0.5*( eta(k) + eta(k+1) )        ! Z of middle of layer (in m)
      t(i,j,k) = sst_ref + dtdz*min(0., zc+mld)
      s(i,j,k) = s_ref
    enddo ! k
  enddo ; enddo

Here is the caller graph for this function:

scm_idealized_hurricane_wind_forcing()

subroutine, public scm_idealized_hurricane::scm_idealized_hurricane_wind_forcing	(	type(surface), intent(in)	state,
		type(forcing), intent(inout)	fluxes,
		type(time_type), intent(in)	day,
		type(ocean_grid_type), intent(inout)	G,
		type(scm_idealized_hurricane_cs), pointer	CS
	)

Parameters

[in]	state	Surface state structure
[in,out]	fluxes	Surface fluxes structure
[in]	day	Time in days
[in,out]	g	Grid structure
	cs	Container for SCM parameters

Definition at line 159 of file SCM_idealized_hurricane.F90.

References mom_forcing_type::allocate_forcing_type().

  type(surface),                    intent(in)    :: state  !< Surface state structure
  type(forcing),                    intent(inout) :: fluxes !< Surface fluxes structure
  type(time_type),                  intent(in)    :: day    !< Time in days
  type(ocean_grid_type),            intent(inout) :: g      !< Grid structure
  type(scm_idealized_hurricane_cs), pointer       :: cs     !< Container for SCM parameters
  ! Local variables
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  real :: pie, deg2rad
  real :: u10, a, b, c, r, f,du10,rkm ! For wind profile expression
  real :: xx, t0 !for location
  real :: dp, rb
  real :: cd ! Air-sea drag coefficient
  real :: uocn, vocn ! Surface ocean velocity components
  real :: du, dv ! Air-sea differential motion
  !Wind angle variables
  real :: alph,rstr, a0, a1, p1, adir, transdir, v_ts, u_ts
  logical :: br_bench
  ! Bounds for loops and memory allocation
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  ! Allocate the forcing arrays, if necessary.
  call allocate_forcing_type(g, fluxes, stress=.true., ustar=.true.)

  pie = 4.0*atan(1.0) ; deg2rad = pie/180.

  !/ BR
  ! Implementing Holland (1980) parameteric wind profile
  !------------------------------------------------------|
  br_bench = .true.   !true if comparing to LES runs     |
  t0 = 129600.        !TC 'eye' crosses (0,0) at 36 hours|
  transdir = pie      !translation direction (-x)        |
  !------------------------------------------------------|
  dp = cs%p_n - cs%p_c
  c = cs%U_max / sqrt( dp )
  b = c**2 * cs%rho_a * exp(1.0)
  if (br_bench) then
     ! rho_a reset to value used in generated wind for benchmark test
     b = c**2 * 1.2 * exp(1.0)
  endif
  a = (cs%r_max/1000.)**b
  f =g%CoriolisBu(is,js) ! f=f(x,y) but in the SCM is constant
  if (br_bench) then
     ! f reset to value used in generated wind for benchmark test
     f = 5.5659e-05
  endif
  !/ BR
  ! Calculate x position as a function of time.
  xx = ( t0 - time_type_to_real(day)) * cs%tran_speed * cos(transdir)
  r = sqrt(xx**2.+cs%YY**2.)
  !/ BR
  ! rkm - r converted to km for Holland prof.
  !       used in km due to error, correct implementation should
  !       not need rkm, but to match winds w/ experiment this must
  !       be maintained.  Causes winds far from storm center to be a
  !       couple of m/s higher than the correct Holland prof.
  if (br_bench) then
     rkm = r/1000.
     rb = (rkm)**b
  else
     ! if not comparing to benchmark, then use correct Holland prof.
     rkm = r
     rb = r**b
  endif
  !/ BR
  ! Calculate U10 in the interior (inside of 10x radius of maximum wind),
  ! while adjusting U10 to 0 outside of 12x radius of maximum wind.
  ! Note that rho_a is set to 1.2 following generated wind for experiment
  if (r/cs%r_max.gt.0.001 .AND. r/cs%r_max.lt.10.) then
     u10 = sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f
  elseif (r/cs%r_max.gt.10. .AND. r/cs%r_max.lt.12.) then
     r=cs%r_max*10.
     if (br_bench) then
        rkm = r/1000.
        rb=rkm**b
     else
        rkm = r
        rb = r**b
     endif
     u10 = ( sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f) &
           * (12. - r/cs%r_max)/2.
  else
     u10 = 0.
  end if
  adir = atan2(cs%YY,xx)

  !/ BR
  ! Wind angle model following Zhang and Ulhorn (2012)
  ! ALPH is inflow angle positive outward.
  rstr = min(10.,r / cs%r_max)
  a0 = -0.9*rstr -0.09*cs%U_max -14.33
  a1 = -a0 *(0.04*rstr +0.05*cs%tran_speed+0.14)
  p1 = (6.88*rstr -9.60*cs%tran_speed+85.31)*pie/180.
  alph = a0 - a1*cos( (transdir - adir ) - p1)
  if (r/cs%r_max.gt.10. .AND. r/cs%r_max.lt.12.) then
     alph = alph* (12. - r/cs%r_max)/2.
  elseif (r/cs%r_max.gt.12.) then
     alph = 0.0
  endif
  alph = alph * deg2rad

  !/BR
  ! Prepare for wind calculation
  ! X_TS is component of translation speed added to wind vector
  ! due to background steering wind.
  u_ts = cs%tran_speed/2.*cos(transdir)
  v_ts = cs%tran_speed/2.*sin(transdir)

  ! Set the surface wind stresses, in units of Pa. A positive taux
  ! accelerates the ocean to the (pseudo-)east.
  !   The i-loop extends to is-1 so that taux can be used later in the
  ! calculation of ustar - otherwise the lower bound would be Isq.
  do j=js,je ; do i=is-1,ieq
    !/BR
    ! Turn off surface current for stress calculation to be
    ! consistent with test case.
    uocn = 0.!state%u(I,j)
    vocn = 0.!0.25*( (state%v(i,J) + state%v(i+1,J-1)) &
             !    +(state%v(i+1,J) + state%v(i,J-1)) )
    !/BR
    ! Wind vector calculated from location/direction (sin/cos flipped b/c
    ! cyclonic wind is 90 deg. phase shifted from position angle).
    du = u10*sin(adir-pie-alph) - uocn + u_ts
    dv = u10*cos(adir-alph) - vocn + v_ts
    !/----------------------------------------------------|
    !BR
    !  Add a simple drag coefficient as a function of U10 |
    !/----------------------------------------------------|
    du10=sqrt(du**2+dv**2)
    if (du10.LT.11.) then
       cd = 1.2e-3
    elseif (du10.LT.20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    fluxes%taux(i,j) = cs%rho_a * g%mask2dCu(i,j) * cd*sqrt(du**2+dv**2)*du
  enddo ; enddo
  !/BR
  ! See notes above
  do j=js-1,jeq ; do i=is,ie
    uocn = 0.!0.25*( (state%u(I,j) + state%u(I-1,j+1)) &
             !    +(state%u(I-1,j) + state%u(I,j+1)) )
    vocn = 0.!state%v(i,J)
    du = u10*sin(adir-pie-alph) - uocn + u_ts
    dv = u10*cos(adir-alph) - vocn + v_ts
    du10=sqrt(du**2+dv**2)
    if (du10.LT.11.) then
       cd = 1.2e-3
    elseif (du10.LT.20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    fluxes%tauy(i,j) = cs%rho_a * g%mask2dCv(i,j) * cd*du10*dv
  enddo ; enddo
  ! Set the surface friction velocity, in units of m s-1. ustar is always positive.
  do j=js,je ; do i=is,ie
    !  This expression can be changed if desired, but need not be.
    fluxes%ustar(i,j) = g%mask2dT(i,j) * sqrt(cs%gust_const/cs%Rho0 + &
       sqrt(0.5*(fluxes%taux(i-1,j)**2 + fluxes%taux(i,j)**2) + &
            0.5*(fluxes%tauy(i,j-1)**2 + fluxes%tauy(i,j)**2))/cs%Rho0)
  enddo ; enddo

Here is the call graph for this function:

scm_idealized_hurricane_wind_init()

subroutine, public scm_idealized_hurricane::scm_idealized_hurricane_wind_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		type(scm_idealized_hurricane_cs), pointer	CS
	)

Initializes wind profile for the SCM idealized hurricane example.

Parameters

[in]	g	Grid structure
[in]	param_file	Input parameter structure
	cs	Parameter container

Definition at line 97 of file SCM_idealized_hurricane.F90.

References mdl, and mom_error_handler::mom_error().

  type(time_type),              intent(in) :: time       !< Time
  type(ocean_grid_type),        intent(in) :: g          !< Grid structure
  type(param_file_type),        intent(in) :: param_file !< Input parameter structure
  type(scm_idealized_hurricane_cs), pointer    :: cs         !< Parameter container

! This include declares and sets the variable "version".
#include "version_variable.h"

  if (associated(cs)) then
    call mom_error(fatal, "SCM_idealized_hurricane_wind_init called with an associated "// &
                          "control structure.")
    return
  endif
  allocate(cs)

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SCM_RHO_AIR", cs%rho_a,            &
                 "Air density "//                                     &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='kg/m3', default=1.2)
  call get_param(param_file, mdl, "SCM_AMBIENT_PRESSURE", cs%p_n,     &
                 "Ambient pressure "//                                &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='Pa', default=101200.)
  call get_param(param_file, mdl, "SCM_CENTRAL_PRESSURE", cs%p_c,     &
                 "Central pressure "//                                &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='Pa', default=96800.)
  call get_param(param_file, mdl, "SCM_RADIUS_MAX_WINDS", cs%r_max,   &
                 "Radius of maximum winds "//                         &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m', default=50.e3)
  call get_param(param_file, mdl, "SCM_MAX_WIND_SPEED", cs%U_max,     &
                 "Maximum wind speed "//                              &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m/s', default=65.)
  call get_param(param_file, mdl, "SCM_YY", cs%YY,     &
                 "Y distance of station "//                           &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m', default=50.e3)
  call get_param(param_file, mdl, "SCM_TRAN_SPEED", cs%TRAN_SPEED,     &
                 "Translation speed of hurricane"//                   &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m/s', default=5.0)
  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to \n"//&
                 "calculate accelerations and the mass for conservation \n"//&
                 "properties, or with BOUSSINSEQ false to convert some \n"//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0)
  ! The following parameter is a model run-time parameter which is used
  ! and logged elsewhere and so should not be logged here. The default
  ! value should be consistent with the rest of the model.
  call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
                 "The background gustiness in the winds.", units="Pa", &
                 default=0.00, do_not_log=.true.)

Here is the call graph for this function:

Variable Documentation

mdl

character(len=40) scm_idealized_hurricane::mdl = "SCM_idealized_hurricane"

private

Definition at line 42 of file SCM_idealized_hurricane.F90.

Referenced by scm_idealized_hurricane_ts_init(), and scm_idealized_hurricane_wind_init().

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