22 implicit none ;
private
24 #include <MOM_memory.h>
42 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
45 real,
intent(in) :: max_depth
52 # include "version_variable.h"
53 character(len=40) :: mdl =
"DOME_initialize_topography"
54 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
55 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
56 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
58 call mom_mesg(
" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5)
60 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
63 call get_param(param_file, mdl,
"MINIMUM_DEPTH", min_depth, &
64 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=m_to_z)
66 do j=js,je ;
do i=is,ie
67 if (g%geoLatT(i,j) < 600.0)
then
68 if (g%geoLatT(i,j) < 300.0)
then
71 d(i,j) = max_depth - 10.0*m_to_z * (g%geoLatT(i,j)-300.0)
74 if ((g%geoLonT(i,j) > 1000.0) .AND. (g%geoLonT(i,j) < 1100.0))
then
77 d(i,j) = 0.5*min_depth
81 if (d(i,j) > max_depth) d(i,j) = max_depth
82 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
93 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
97 logical,
optional,
intent(in) :: just_read_params
100 real :: e0(szk_(gv)+1)
102 real :: eta1d(szk_(gv)+1)
105 character(len=40) :: mdl =
"DOME_initialize_thickness"
106 integer :: i, j, k, is, ie, js, je, nz
108 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
110 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
112 if (just_read)
return
114 call mom_mesg(
" DOME_initialization.F90, DOME_initialize_thickness: setting thickness", 5)
118 e0(k) = -g%max_depth * (real(k-1)-0.5)/real(nz-1)
121 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
127 eta1d(nz+1) = -g%bathyT(i,j)
130 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
131 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
132 h(i,j,k) = gv%Angstrom_H
134 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
160 real :: eta(szi_(g),szj_(g),szk_(g)+1)
161 real :: temp(szi_(g),szj_(g),szk_(g))
162 real :: idamp(szi_(g),szj_(g))
166 real :: damp, e_dense, damp_new
167 character(len=40) :: mdl =
"DOME_initialize_sponges"
168 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
170 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
171 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
173 eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
181 call get_param(pf, mdl,
"MINIMUM_DEPTH", min_depth, &
182 "The minimum depth of the ocean.", units=
"m", default=0.0, scale=us%m_to_Z)
185 do k=2,nz ; h0(k) = -(real(k-1)-0.5)*g%max_depth / real(nz-1) ;
enddo
186 do i=is,ie;
do j=js,je
187 if (g%geoLonT(i,j) < 100.0)
then ; damp = 10.0
188 elseif (g%geoLonT(i,j) < 200.0)
then
189 damp = 10.0*(200.0-g%geoLonT(i,j))/100.0
193 if (g%geoLonT(i,j) > 1400.0)
then ; damp_new = 10.0
194 elseif (g%geoLonT(i,j) > 1300.0)
then
195 damp_new = 10.0*(g%geoLonT(i,j)-1300.0)/100.0
196 else ; damp_new = 0.0
199 if (damp <= damp_new) damp=damp_new
206 e_dense = -g%bathyT(i,j)
207 if (e_dense >= h0(k))
then ; eta(i,j,k) = e_dense
208 else ; eta(i,j,k) = h0(k) ;
endif
209 if (eta(i,j,k) < gv%Angstrom_Z*(nz-k+1) - g%bathyT(i,j)) &
210 eta(i,j,k) = gv%Angstrom_Z*(nz-k+1) - g%bathyT(i,j)
212 eta(i,j,nz+1) = -g%bathyT(i,j)
214 if (g%bathyT(i,j) > min_depth)
then
215 idamp(i,j) = damp/86400.0
216 else ; idamp(i,j) = 0.0 ;
endif
231 if (
associated(tv%T) )
then
232 call mom_error(fatal,
"DOME_initialize_sponges is not set up for use with"//&
233 " a temperatures defined.")
261 real :: t0(szk_(g)), s0(szk_(g))
262 real :: pres(szk_(g))
263 real :: drho_dt(szk_(g))
264 real :: drho_ds(szk_(g))
265 real :: rho_guess(szk_(g))
267 real :: tr_0, y1, y2, tr_k, rst, rsb, rc, v_k, lon_im1
275 character(len=40) :: mdl =
"DOME_set_OBC_data"
276 character(len=32) :: name
277 integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntr
278 integer :: isdb, iedb, jsdb, jedb
282 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
283 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
284 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
287 d_edge = 300.0*us%m_to_Z
291 if (.not.
associated(obc))
return
293 g_prime_tot = (gv%g_Earth / gv%Rho0) * 2.0*us%kg_m3_to_R
294 def_rad = us%L_to_m*sqrt(d_edge*g_prime_tot) / (1.0e-4*us%T_to_s * 1000.0)
295 tr_0 = (-d_edge*sqrt(d_edge*g_prime_tot)*0.5e3*us%m_to_L*def_rad) * gv%Z_to_H
297 if (obc%number_of_segments /= 1)
then
298 call mom_error(warning,
'Error in DOME OBC segment setup', .true.)
306 if (.not.
associated(obc%tracer_x_reservoirs_used))
then
307 allocate(obc%tracer_x_reservoirs_used(ntr))
308 allocate(obc%tracer_y_reservoirs_used(ntr))
309 obc%tracer_x_reservoirs_used(:) = .false.
310 obc%tracer_y_reservoirs_used(:) = .false.
311 obc%tracer_y_reservoirs_used(1) = .true.
314 segment => obc%segment(1)
315 if (.not. segment%on_pe)
return
317 allocate(segment%field(ntr))
321 if (k>1) rst = -1.0 + (real(k-1)-0.5)/real(nz-1)
324 if (k<nz) rsb = -1.0 + (real(k-1)+0.5)/real(nz-1)
325 rc = -1.0 + real(k-1)/real(nz-1)
328 y1 = (2.0*ri_trans*rst + ri_trans + 2.0)/(2.0 - ri_trans)
329 y2 = (2.0*ri_trans*rsb + ri_trans + 2.0)/(2.0 - ri_trans)
330 tr_k = tr_0 * (2.0/(ri_trans*(2.0-ri_trans))) * &
331 ((log(y1)+1.0)/y1 - (log(y2)+1.0)/y2)
332 v_k = -sqrt(d_edge*g_prime_tot)*log((2.0 + ri_trans*(1.0 + 2.0*rc)) / &
334 if (k == nz) tr_k = tr_k + tr_0 * (2.0/(ri_trans*(2.0+ri_trans))) * &
335 log((2.0+ri_trans)/(2.0-ri_trans))
337 isd = segment%HI%isd ; ied = segment%HI%ied
338 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
339 do j=jsdb,jedb ;
do i=isd,ied
340 lon_im1 = 2.0*g%geoLonCv(i,j) - g%geoLonBu(i,j)
341 segment%normal_trans(i,j,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/def_rad) -&
342 exp(-2.0*(g%geoLonBu(i,j) - 1000.0)/def_rad))
343 segment%normal_vel(i,j,k) = v_k * exp(-2.0*(g%geoLonCv(i,j) - 1000.0)/def_rad)
349 if (
associated(tv%S))
then
355 if (
associated(tv%T))
then
359 pres(:) = tv%P_Ref ; s0(:) = 35.0 ; t0(1) = 25.0
360 call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),tv%eqn_of_state, scale=us%kg_m3_to_R)
363 do k=1,nz ; t0(k) = t0(1) + (gv%Rlay(k)-rho_guess(1)) / drho_dt(1) ;
enddo
365 call calculate_density(t0,s0,pres,rho_guess,1,nz,tv%eqn_of_state, scale=us%kg_m3_to_R)
367 do k=1,nz ; t0(k) = t0(k) + (gv%Rlay(k)-rho_guess(k)) / drho_dt(k) ;
enddo
371 allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
372 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied
373 segment%field(1)%buffer_src(i,j,k) = t0(k)
374 enddo ;
enddo ;
enddo
383 allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
384 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%isd,segment%HI%ied
385 if (k < nz/2)
then ; segment%field(1)%buffer_src(i,j,k) = 0.0
386 else ; segment%field(1)%buffer_src(i,j,k) = 1.0 ;
endif
387 enddo ;
enddo ;
enddo
391 obc%segment(1), obc_array=.true.)
396 if (m < 10)
then ;
write(name,
'("tr_D",I1.1)') m
397 else ;
write(name,
'("tr_D",I2.2)') m ;
endif
400 obc%segment(1), obc_scalar=0.0)