18 implicit none ;
private
20 #include <MOM_memory.h>
36 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
39 real,
intent(in) :: max_depth
50 # include "version_variable.h"
51 character(len=40) :: mdl =
"benchmark_initialize_topography"
52 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
53 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
54 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
56 call mom_mesg(
" benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5)
58 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
61 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
62 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=m_to_z)
68 do j=js,je ;
do i=is,ie
69 x = (g%geoLonT(i,j)-g%west_lon) / g%len_lon
70 y = (g%geoLatT(i,j)-g%south_lat) / g%len_lat
72 d(i,j) = -d0 * ( y*(1.0 + 0.6*cos(4.0*pi*x)) &
74 + 0.05*cos(10.0*pi*x) - 0.7 )
75 if (d(i,j) > max_depth) d(i,j) = max_depth
76 if (d(i,j) < min_depth) d(i,j) = 0.
86 P_ref, just_read_params)
90 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
94 type(
eos_type),
pointer :: eqn_of_state
96 real,
intent(in) :: p_ref
98 logical,
optional,
intent(in) :: just_read_params
101 real :: e0(szk_(gv)+1)
103 real :: e_pert(szk_(gv)+1)
105 real :: eta1d(szk_(gv)+1)
110 real :: thermocline_scale
111 real,
dimension(SZK_(GV)) :: &
125 # include "version_variable.h"
126 character(len=40) :: mdl =
"benchmark_initialize_thickness"
127 integer :: i, j, k, k1, is, ie, js, je, nz, itt
129 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
131 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
132 if (.not.just_read)
call log_version(param_file, mdl, version,
"")
133 call get_param(param_file, mdl,
"BENCHMARK_ML_DEPTH_IC", ml_depth, &
134 "Initial mixed layer depth in the benchmark test case.", &
135 units=
'm', default=50.0, scale=us%m_to_Z, do_not_log=just_read)
136 call get_param(param_file, mdl,
"BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, &
137 "Initial thermocline depth scale in the benchmark test case.", &
138 default=500.0, units=
"m", scale=us%m_to_Z, do_not_log=just_read)
140 if (just_read)
return
142 call mom_mesg(
" benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5)
144 k1 = gv%nk_rho_varies + 1
151 pres(k) = p_ref ; s0(k) = 35.0
154 call calculate_density(t0(k1), s0(k1), pres(k1), rho_guess(k1), eqn_of_state, scale=us%kg_m3_to_R)
159 t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
164 call calculate_density(t0, s0, pres, rho_guess, 1, nz, eqn_of_state, scale=us%kg_m3_to_R)
167 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
172 i_ts = 1.0 / thermocline_scale
173 i_md = 1.0 / g%max_depth
174 do j=js,je ;
do i=is,ie
175 sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
176 cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
178 do k=1,nz ; e_pert(k) = 0.0 ;
enddo
184 eta1d(nz+1) = -g%bathyT(i,j)
187 t_int = 0.5*(t0(k) + t0(k-1))
188 t_frac = (t_int - t0(nz)) / (sst - t0(nz))
192 err = a_exp * exp(z*i_ts) + (1.0 - a_exp) * (z*i_md + 1.0) - t_frac
193 derr_dz = a_exp * i_ts * exp(z*i_ts) + (1.0 - a_exp) * i_md
194 z = z - err / derr_dz
199 eta1d(k) = e0(k) + e_pert(k)
201 if (eta1d(k) > -ml_depth) eta1d(k) = -ml_depth
203 if (eta1d(k) < eta1d(k+1) + gv%Angstrom_Z) &
204 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
206 h(i,j,k) = max(gv%Z_to_H * (eta1d(k) - eta1d(k+1)), gv%Angstrom_H)
208 h(i,j,1) = max(gv%Z_to_H * (0.0 - eta1d(2)), gv%Angstrom_H)
216 eqn_of_state, P_Ref, just_read_params)
219 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
221 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
227 type(
eos_type),
pointer :: eqn_of_state
229 real,
intent(in) :: p_ref
231 logical,
optional,
intent(in) :: just_read_params
234 real :: t0(szk_(g)), s0(szk_(g))
235 real :: pres(szk_(g))
236 real :: drho_dt(szk_(g))
237 real :: drho_ds(szk_(g))
238 real :: rho_guess(szk_(g))
243 character(len=40) :: mdl =
"benchmark_init_temperature_salinity"
244 integer :: i, j, k, k1, is, ie, js, je, nz, itt
246 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
248 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
250 if (just_read)
return
252 k1 = gv%nk_rho_varies + 1
255 pres(k) = p_ref ; s0(k) = 35.0
259 call calculate_density(t0(k1),s0(k1),pres(k1),rho_guess(k1),eqn_of_state, scale=us%kg_m3_to_R)
264 t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
269 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
272 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
276 do k=1,nz ;
do i=is,ie ;
do j=js,je
279 enddo ;
enddo ;
enddo
281 do i=is,ie ;
do j=js,je
282 sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
283 cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))