18 implicit none ;
private
20 #include <MOM_memory.h>
23 character(len=40) ::
mdl =
"Rossby_front_2d_initialization"
25 #include "version_variable.h"
43 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
47 logical,
optional,
intent(in) :: just_read_params
50 integer :: i, j, k, is, ie, js, je, nz
51 real :: tz, dml, eta, stretch, h0
52 real :: min_thickness, t_range
55 character(len=40) :: verticalcoordinate
57 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
59 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
62 call mom_mesg(
"Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness")
66 call get_param(param_file,
mdl,
"MIN_THICKNESS", min_thickness, &
67 'Minimum layer thickness',units=
'm',default=1.e-3, do_not_log=just_read)
68 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
69 default=default_coordinate_mode, do_not_log=just_read)
70 call get_param(param_file,
mdl,
"T_RANGE", t_range,
'Initial temperature range', &
71 units=
'C', default=0.0, do_not_log=just_read)
72 call get_param(param_file,
mdl,
"DRHO_DT", drho_dt, default=-0.2, scale=us%kg_m3_to_R, do_not_log=.true.)
76 tz = t_range / g%max_depth
78 select case ( coordinatemode(verticalcoordinate) )
81 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
82 dml =
hml( g, g%geoLatT(i,j) )
83 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
84 stretch = ( ( g%max_depth + eta ) / g%max_depth )
85 h0 = ( g%max_depth / real(nz) ) * stretch
87 h(i,j,k) = h0 * gv%Z_to_H
92 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
93 dml =
hml( g, g%geoLatT(i,j) )
94 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
95 stretch = ( ( g%max_depth + eta ) / g%max_depth )
96 h0 = ( g%max_depth / real(nz) ) * stretch
98 h(i,j,k) = h0 * gv%Z_to_H
103 call mom_error(fatal,
"Rossby_front_initialize: "// &
104 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
113 param_file, eqn_of_state, just_read_params)
116 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
117 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
118 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
120 type(
eos_type),
pointer :: eqn_of_state
121 logical,
optional,
intent(in) :: just_read_params
124 integer :: i, j, k, is, ie, js, je, nz
127 real :: y, zc, zi, dtdz
129 character(len=40) :: verticalcoordinate
132 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
134 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
136 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
137 default=default_coordinate_mode, do_not_log=just_read)
138 call get_param(param_file,
mdl,
"S_REF",s_ref,
'Reference salinity', units=
'1e-3', &
139 fail_if_missing=.not.just_read, do_not_log=just_read)
140 call get_param(param_file,
mdl,
"T_REF",t_ref,
'Reference temperature',units=
'C',&
141 fail_if_missing=.not.just_read, do_not_log=just_read)
142 call get_param(param_file,
mdl,
"T_RANGE",t_range,
'Initial temperature range',&
143 units=
'C', default=0.0, do_not_log=just_read)
145 if (just_read)
return
149 dtdz = t_range / g%max_depth
151 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
155 zc = gv%H_to_Z * (zi - 0.5*h(i,j,k))
156 zc = min( zc, -
hml(g, g%geoLatT(i,j)) )
157 t(i,j,k) = t_ref + dtdz * zc
168 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
170 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
172 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)), &
177 logical,
optional,
intent(in) :: just_read_params
184 real :: dml, zi, zc, zm
188 integer :: i, j, k, is, ie, js, je, nz
190 character(len=40) :: verticalcoordinate
192 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
194 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
196 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
197 default=default_coordinate_mode, do_not_log=just_read)
198 call get_param(param_file,
mdl,
"T_RANGE", t_range,
'Initial temperature range', &
199 units=
'C', default=0.0, do_not_log=just_read)
200 call get_param(param_file,
mdl,
"DRHO_DT", drho_dt, default=-0.2, scale=us%kg_m3_to_R, do_not_log=.true.)
202 if (just_read)
return
207 do j = g%jsc,g%jec ;
do i = g%isc-1,g%iec+1
208 f = 0.5* (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
209 dudt = 0.0 ;
if (abs(f) > 0.0) &
210 dudt = ( gv%g_Earth*drho_dt ) / ( f * gv%Rho0 )
211 dml =
hml( g, g%geoLatT(i,j) )
212 ty = us%L_to_m*
dtdy( g, t_range, g%geoLatT(i,j) )
215 hatu = 0.5*(h(i,j,k)+h(i+1,j,k)) * gv%H_to_Z
218 zm = max( zc + dml, 0. )
219 u(i,j,k) = dudt * ty * zm
228 real function ypseudo( G, lat )
230 real,
intent(in) :: lat
235 ypseudo = ( ( lat - g%south_lat ) / g%len_lat ) - 0.5
242 real function hml( G, lat )
244 real,
intent(in) :: lat
246 real :: dhml, hmlmean
255 real function dtdy( G, dT, lat )
257 real,
intent(in) :: dt
258 real,
intent(in) :: lat
260 real :: pi, dhml, dhdy
266 dtdy = -( dt / g%max_depth ) * dhdy