14 use regrid_consts,
only : coordinatemode, default_coordinate_mode
18 implicit none ;
private
20 #include <MOM_memory.h>
23 character(len=40) :: mdl =
"Rossby_front_2d_initialization"
25 #include "version_variable.h"
27 public rossby_front_initialize_thickness
28 public rossby_front_initialize_temperature_salinity
29 public rossby_front_initialize_velocity
32 real,
parameter :: frontfractionalwidth = 0.5
33 real,
parameter :: hmlmin = 0.25
34 real,
parameter :: hmlmax = 0.75
39 subroutine rossby_front_initialize_thickness(h, G, GV, param_file, just_read_params)
42 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
46 logical,
optional,
intent(in) :: just_read_params
49 integer :: i, j, k, is, ie, js, je, nz
50 real :: tz, dml, eta, stretch, h0
51 real :: min_thickness, t_range, drho_dt
53 character(len=40) :: verticalcoordinate
55 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
57 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
60 call mom_mesg(
"Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness")
62 if (.not.just_read)
call log_version(param_file, mdl, version,
"")
64 call get_param(param_file, mdl,
"MIN_THICKNESS", min_thickness, &
65 'Minimum layer thickness',units=
'm',default=1.e-3, do_not_log=just_read)
66 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
67 default=default_coordinate_mode, do_not_log=just_read)
68 call get_param(param_file, mdl,
"T_RANGE", t_range,
'Initial temperature range', &
69 units=
'C', default=0.0, do_not_log=just_read)
70 call get_param(param_file, mdl,
"DRHO_DT", drho_dt, default=-0.2, do_not_log=.true.)
74 tz = t_range / g%max_depth
76 select case ( coordinatemode(verticalcoordinate) )
78 case (regridding_layer, regridding_rho)
79 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
80 dml = hml( g, g%geoLatT(i,j) )
81 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
82 stretch = ( ( g%max_depth + eta ) / g%max_depth )
83 h0 = ( g%max_depth / real(nz) ) * stretch
85 h(i,j,k) = h0 * gv%Z_to_H
89 case (regridding_zstar, regridding_sigma)
90 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
91 dml = hml( g, g%geoLatT(i,j) )
92 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
93 stretch = ( ( g%max_depth + eta ) / g%max_depth )
94 h0 = ( g%max_depth / real(nz) ) * stretch
96 h(i,j,k) = h0 * gv%Z_to_H
101 call mom_error(fatal,
"Rossby_front_initialize: "// &
102 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
106 end subroutine rossby_front_initialize_thickness
110 subroutine rossby_front_initialize_temperature_salinity(T, S, h, G, GV, &
111 param_file, eqn_of_state, just_read_params)
114 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
115 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
116 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
118 type(
eos_type),
pointer :: eqn_of_state
119 logical,
optional,
intent(in) :: just_read_params
122 integer :: i, j, k, is, ie, js, je, nz
125 real :: y, zc, zi, dtdz
127 character(len=40) :: verticalcoordinate
130 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
132 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
134 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
135 default=default_coordinate_mode, do_not_log=just_read)
136 call get_param(param_file, mdl,
"S_REF",s_ref,
'Reference salinity', units=
'1e-3', &
137 fail_if_missing=.not.just_read, do_not_log=just_read)
138 call get_param(param_file, mdl,
"T_REF",t_ref,
'Reference temperature',units=
'C',&
139 fail_if_missing=.not.just_read, do_not_log=just_read)
140 call get_param(param_file, mdl,
"T_RANGE",t_range,
'Initial temperature range',&
141 units=
'C', default=0.0, do_not_log=just_read)
143 if (just_read)
return
147 dtdz = t_range / g%max_depth
149 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
153 zc = gv%H_to_Z * (zi - 0.5*h(i,j,k))
154 zc = min( zc, -hml(g, g%geoLatT(i,j)) )
155 t(i,j,k) = t_ref + dtdz * zc
159 end subroutine rossby_front_initialize_temperature_salinity
163 subroutine rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just_read_params)
167 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
169 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
171 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)), &
175 logical,
optional,
intent(in) :: just_read_params
182 real :: dml, zi, zc, zm
186 integer :: i, j, k, is, ie, js, je, nz
188 character(len=40) :: verticalcoordinate
190 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
192 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
194 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
195 default=default_coordinate_mode, do_not_log=just_read)
196 call get_param(param_file, mdl,
"T_RANGE", t_range,
'Initial temperature range', &
197 units=
'C', default=0.0, do_not_log=just_read)
198 call get_param(param_file, mdl,
"DRHO_DT", drho_dt, default=-0.2, do_not_log=.true.)
200 if (just_read)
return
205 do j = g%jsc,g%jec ;
do i = g%isc-1,g%iec+1
206 f = 0.5* (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
207 dudt = 0.0 ;
if (abs(f) > 0.0) &
208 dudt = ( gv%g_Earth*drho_dt ) / ( f * gv%Rho0 )
209 dml = hml( g, g%geoLatT(i,j) )
210 ty = us%L_to_m*dtdy( g, t_range, g%geoLatT(i,j) )
213 hatu = 0.5*(h(i,j,k)+h(i+1,j,k)) * gv%H_to_Z
216 zm = max( zc + dml, 0. )
217 u(i,j,k) = us%L_T_to_m_s * dudt * ty * zm
221 end subroutine rossby_front_initialize_velocity
226 real function ypseudo( G, lat )
228 real,
intent(in) :: lat
233 ypseudo = ( ( lat - g%south_lat ) / g%len_lat ) - 0.5
234 ypseudo = pi * max(-0.5, min(0.5, ypseudo / frontfractionalwidth))
240 real function hml( G, lat )
242 real,
intent(in) :: lat
244 real :: dhml, hmlmean
246 dhml = 0.5 * ( hmlmax - hmlmin ) * g%max_depth
247 hmlmean = 0.5 * ( hmlmin + hmlmax ) * g%max_depth
248 hml = hmlmean + dhml * sin( ypseudo(g, lat) )
253 real function dtdy( G, dT, lat )
255 real,
intent(in) :: dt
256 real,
intent(in) :: lat
258 real :: pi, dhml, dhdy
262 dhml = 0.5 * ( hmlmax - hmlmin ) * g%max_depth
263 dhdy = dhml * ( pi / ( frontfractionalwidth * g%len_lat * km ) ) * cos( ypseudo(g, lat) )
264 dtdy = -( dt / g%max_depth ) * dhdy