22 implicit none ;
private
24 #include <MOM_memory.h>
26 public dome_initialize_topography
27 public dome_initialize_thickness
28 public dome_initialize_sponges
29 public dome_set_obc_data
40 subroutine dome_initialize_topography(D, G, param_file, max_depth, US)
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
85 end subroutine dome_initialize_topography
90 subroutine dome_initialize_thickness(h, G, GV, param_file, just_read_params)
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))
139 end subroutine dome_initialize_thickness
148 subroutine dome_initialize_sponges(G, GV, US, tv, PF, CSp)
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
221 call initialize_sponge(idamp, eta, g, pf, csp, gv)
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.")
235 call set_up_sponge_field(temp, tv%T, g, nz, csp)
237 call set_up_sponge_field(temp, tv%S, g, nz, csp)
240 end subroutine dome_initialize_sponges
244 subroutine dome_set_obc_data(OBC, tv, G, GV, US, param_file, tr_Reg)
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
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%s_to_T*us%L_to_m*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.)
301 segment => obc%segment(1)
302 if (.not. segment%on_pe)
return
305 allocate(segment%field(ntr))
309 if (k>1) rst = -1.0 + (real(k-1)-0.5)/real(nz-1)
312 if (k<nz) rsb = -1.0 + (real(k-1)+0.5)/real(nz-1)
313 rc = -1.0 + real(k-1)/real(nz-1)
316 y1 = (2.0*ri_trans*rst + ri_trans + 2.0)/(2.0 - ri_trans)
317 y2 = (2.0*ri_trans*rsb + ri_trans + 2.0)/(2.0 - ri_trans)
318 tr_k = tr_0 * (2.0/(ri_trans*(2.0-ri_trans))) * &
319 ((log(y1)+1.0)/y1 - (log(y2)+1.0)/y2)
320 v_k = -us%L_T_to_m_s*sqrt(d_edge*g_prime_tot)*log((2.0 + ri_trans*(1.0 + 2.0*rc)) / &
322 if (k == nz) tr_k = tr_k + tr_0 * (2.0/(ri_trans*(2.0+ri_trans))) * &
323 log((2.0+ri_trans)/(2.0-ri_trans))
325 isd = segment%HI%isd ; ied = segment%HI%ied
326 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
327 do j=jsdb,jedb ;
do i=isd,ied
328 lon_im1 = 2.0*g%geoLonCv(i,j) - g%geoLonBu(i,j)
329 segment%normal_trans(i,j,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/def_rad) -&
330 exp(-2.0*(g%geoLonBu(i,j) - 1000.0)/def_rad))
331 segment%normal_vel(i,j,k) = v_k * exp(-2.0*(g%geoLonCv(i,j) - 1000.0)/def_rad)
337 if (
associated(tv%S))
then
340 call tracer_name_lookup(tr_reg, tr_ptr, name)
341 call register_segment_tracer(tr_ptr, param_file, gv, segment, obc_scalar=35.0)
343 if (
associated(tv%T))
then
347 pres(:) = tv%P_Ref ; s0(:) = 35.0 ; t0(1) = 25.0
351 do k=1,nz ; t0(k) = t0(1) + (gv%Rlay(k)-rho_guess(1)) / drho_dt(1) ;
enddo
355 do k=1,nz ; t0(k) = t0(k) + (gv%Rlay(k)-rho_guess(k)) / drho_dt(k) ;
enddo
359 allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
360 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied
361 segment%field(1)%buffer_src(i,j,k) = t0(k)
362 enddo ;
enddo ;
enddo
364 call tracer_name_lookup(tr_reg, tr_ptr, name)
365 call register_segment_tracer(tr_ptr, param_file, gv, segment, obc_array=.true.)
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=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%isd,segment%HI%ied
373 if (k < nz/2)
then ; segment%field(1)%buffer_src(i,j,k) = 0.0
374 else ; segment%field(1)%buffer_src(i,j,k) = 1.0 ;
endif
375 enddo ;
enddo ;
enddo
377 call tracer_name_lookup(tr_reg, tr_ptr, name)
378 call register_segment_tracer(tr_ptr, param_file, gv, &
379 obc%segment(1), obc_array=.true.)
384 if (m < 10)
then ;
write(name,
'("tr_D",I1.1)') m
385 else ;
write(name,
'("tr_D",I2.2)') m ;
endif
386 call tracer_name_lookup(tr_reg, tr_ptr, name)
387 call register_segment_tracer(tr_ptr, param_file, gv, &
388 obc%segment(1), obc_scalar=0.0)
391 end subroutine dome_set_obc_data