18 implicit none ;
private
20 #include <MOM_memory.h>
22 public benchmark_initialize_topography
23 public benchmark_initialize_thickness
24 public benchmark_init_temperature_salinity
34 subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US)
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.
79 end subroutine benchmark_initialize_topography
85 subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, &
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)) :: t0, pres, s0, rho_guess, drho, drho_dt, drho_ds
121 # include "version_variable.h"
122 character(len=40) :: mdl =
"benchmark_initialize_thickness"
123 integer :: i, j, k, k1, is, ie, js, je, nz, itt
125 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
127 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
128 if (.not.just_read)
call log_version(param_file, mdl, version,
"")
129 call get_param(param_file, mdl,
"BENCHMARK_ML_DEPTH_IC", ml_depth, &
130 "Initial mixed layer depth in the benchmark test case.", &
131 units=
'm', default=50.0, scale=us%m_to_Z, do_not_log=just_read)
132 call get_param(param_file, mdl,
"BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, &
133 "Initial thermocline depth scale in the benchmark test case.", &
134 default=500.0, units=
"m", scale=us%m_to_Z, do_not_log=just_read)
136 if (just_read)
return
138 call mom_mesg(
" benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5)
140 k1 = gv%nk_rho_varies + 1
147 pres(k) = p_ref ; s0(k) = 35.0
155 t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
163 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
168 i_ts = 1.0 / thermocline_scale
169 i_md = 1.0 / g%max_depth
170 do j=js,je ;
do i=is,ie
171 sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
172 cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
174 do k=1,nz ; e_pert(k) = 0.0 ;
enddo
180 eta1d(nz+1) = -g%bathyT(i,j)
183 t_int = 0.5*(t0(k) + t0(k-1))
184 t_frac = (t_int - t0(nz)) / (sst - t0(nz))
188 err = a_exp * exp(z*i_ts) + (1.0 - a_exp) * (z*i_md + 1.0) - t_frac
189 derr_dz = a_exp * i_ts * exp(z*i_ts) + (1.0 - a_exp) * i_md
190 z = z - err / derr_dz
195 eta1d(k) = e0(k) + e_pert(k)
197 if (eta1d(k) > -ml_depth) eta1d(k) = -ml_depth
199 if (eta1d(k) < eta1d(k+1) + gv%Angstrom_Z) &
200 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
202 h(i,j,k) = max(gv%Z_to_H * (eta1d(k) - eta1d(k+1)), gv%Angstrom_H)
204 h(i,j,1) = max(gv%Z_to_H * (0.0 - eta1d(2)), gv%Angstrom_H)
208 end subroutine benchmark_initialize_thickness
211 subroutine benchmark_init_temperature_salinity(T, S, G, GV, param_file, &
212 eqn_of_state, P_Ref, just_read_params)
215 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
217 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
222 type(
eos_type),
pointer :: eqn_of_state
224 real,
intent(in) :: p_ref
226 logical,
optional,
intent(in) :: just_read_params
229 real :: t0(szk_(g)), s0(szk_(g))
230 real :: pres(szk_(g))
231 real :: drho_dt(szk_(g))
232 real :: drho_ds(szk_(g))
233 real :: rho_guess(szk_(g))
238 character(len=40) :: mdl =
"benchmark_init_temperature_salinity"
239 integer :: i, j, k, k1, is, ie, js, je, nz, itt
241 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
243 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
245 if (just_read)
return
247 k1 = gv%nk_rho_varies + 1
250 pres(k) = p_ref ; s0(k) = 35.0
259 t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
267 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
271 do k=1,nz ;
do i=is,ie ;
do j=js,je
274 enddo ;
enddo ;
enddo
276 do i=is,ie ;
do j=js,je
277 sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
278 cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
284 end subroutine benchmark_init_temperature_salinity