23 use coupler_types_mod,
only : coupler_type_set_data, ind_csurf
26 implicit none ;
private
28 #include <MOM_memory.h>
34 integer,
parameter ::
ntr = 11
39 logical :: coupled_tracers = .false.
40 character(len=200) :: tracer_ic_file
41 type(time_type),
pointer :: time => null()
43 real,
pointer :: tr(:,:,:,:) => null()
44 real :: land_val(
ntr) = -1.0
46 logical :: tracers_may_reinit
54 integer,
dimension(NTR) :: ind_tr
79 character(len=80) :: name, longname
81 #include "version_variable.h"
82 character(len=40) :: mdl =
"advection_test_tracer"
83 character(len=200) :: inputdir
84 character(len=48) :: flux_units
86 real,
pointer :: tr_ptr(:,:,:) => null()
88 integer :: isd, ied, jsd, jed, nz, m
89 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
91 if (
associated(cs))
then
92 call mom_error(warning,
"register_advection_test_tracer called with an "// &
93 "associated control structure.")
101 call get_param(param_file, mdl,
"ADVECTION_TEST_X_ORIGIN", cs%x_origin, &
102 "The x-coorindate of the center of the test-functions.", default=0.)
103 call get_param(param_file, mdl,
"ADVECTION_TEST_Y_ORIGIN", cs%y_origin, &
104 "The y-coorindate of the center of the test-functions.", default=0.)
105 call get_param(param_file, mdl,
"ADVECTION_TEST_X_WIDTH", cs%x_width, &
106 "The x-width of the test-functions.", default=0.)
107 call get_param(param_file, mdl,
"ADVECTION_TEST_Y_WIDTH", cs%y_width, &
108 "The y-width of the test-functions.", default=0.)
109 call get_param(param_file, mdl,
"ADVECTION_TEST_TRACER_IC_FILE", cs%tracer_IC_file, &
110 "The name of a file from which to read the initial "//&
111 "conditions for the tracers, or blank to initialize "//&
112 "them internally.", default=
" ")
114 if (len_trim(cs%tracer_IC_file) >= 1)
then
115 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
116 cs%tracer_IC_file = trim(slasher(inputdir))//trim(cs%tracer_IC_file)
117 call log_param(param_file, mdl,
"INPUTDIR/ADVECTION_TEST_TRACER_IC_FILE", &
120 call get_param(param_file, mdl,
"SPONGE", cs%use_sponge, &
121 "If true, sponges may be applied anywhere in the domain. "//&
122 "The exact location and properties of those sponges are "//&
123 "specified from MOM_initialization.F90.", default=.false.)
125 call get_param(param_file, mdl,
"TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
126 "If true, tracers may go through the initialization code "//&
127 "if they are not found in the restart files. Otherwise "//&
128 "it is a fatal error if the tracers are not found in the "//&
129 "restart files of a restarted run.", default=.false.)
132 allocate(cs%tr(isd:ied,jsd:jed,nz,
ntr)) ; cs%tr(:,:,:,:) = 0.0
135 if (m < 10)
then ;
write(name,
'("tr",I1.1)') m
136 else ;
write(name,
'("tr",I2.2)') m ;
endif
137 write(longname,
'("Concentration of Tracer ",I2.2)') m
138 cs%tr_desc(m) =
var_desc(name, units=
"kg kg-1", longname=longname, caller=mdl)
139 if (gv%Boussinesq)
then ; flux_units =
"kg kg-1 m3 s-1"
140 else ; flux_units =
"kg s-1" ;
endif
145 tr_ptr => cs%tr(:,:,:,m)
148 name=name, longname=longname, units=
"kg kg-1", &
149 registry_diags=.true., flux_units=flux_units, &
150 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
155 if (cs%coupled_tracers) &
157 flux_type=
' ', implementation=
' ', caller=
"register_advection_test_tracer")
161 cs%restart_CSp => restart_cs
168 logical,
intent(in) :: restart
170 type(time_type),
target,
intent(in) :: day
173 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
175 type(
diag_ctrl),
target,
intent(in) :: diag
185 real,
allocatable :: temp(:,:,:)
186 real,
pointer,
dimension(:,:,:) :: &
187 obc_tr1_u => null(), &
191 character(len=16) :: name
192 character(len=72) :: longname
193 character(len=48) :: units
194 character(len=48) :: flux_units
196 real,
pointer :: tr_ptr(:,:,:) => null()
199 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
200 integer :: isdb, iedb, jsdb, jedb
201 real :: tmpx, tmpy, locx, locy
203 if (.not.
associated(cs))
return
204 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
205 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
206 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
207 h_neglect = gv%H_subroundoff
213 caller=
"initialize_advection_test_tracer")
214 if ((.not.restart) .or. (cs%tracers_may_reinit .and. .not. &
216 do k=1,nz ;
do j=js,je ;
do i=is,ie
218 enddo ;
enddo ;
enddo
220 do j=js,je ;
do i=is,ie
221 if (abs(g%geoLonT(i,j)-cs%x_origin)<0.5*cs%x_width .and. &
222 abs(g%geoLatT(i,j)-cs%y_origin)<0.5*cs%y_width) cs%tr(i,j,k,m) = 1.0
225 do j=js,je ;
do i=is,ie
226 locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
227 locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
228 cs%tr(i,j,k,m) = max(0.0, 1.0-locx)*max(0.0, 1.0-locy)
231 do j=js,je ;
do i=is,ie
232 locx=min(1.0, abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width)*(acos(0.0)*2.)
233 locy=min(1.0, abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width)*(acos(0.0)*2.)
234 cs%tr(i,j,k,m) = (1.0+cos(locx))*(1.0+cos(locy))*0.25
237 do j=js,je ;
do i=is,ie
238 locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
239 locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
240 if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
243 do j=js,je ;
do i=is,ie
244 locx=(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
245 locy=(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
246 if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
247 if (locx>0.0.and.abs(locy)<0.2) cs%tr(i,j,k,m) = 0.0
258 subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
259 evap_CFL_limit, minimum_forcing_depth)
262 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
264 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
266 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
270 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
274 type(
forcing),
intent(in) :: fluxes
276 real,
intent(in) :: dt
280 real,
optional,
intent(in) :: evap_cfl_limit
282 real,
optional,
intent(in) :: minimum_forcing_depth
291 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
293 real :: c1(szi_(g),szk_(g))
294 integer :: i, j, k, is, ie, js, je, nz, m
295 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
297 if (.not.
associated(cs))
return
299 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
301 do k=1,nz ;
do j=js,je ;
do i=is,ie
302 h_work(i,j,k) = h_old(i,j,k)
303 enddo ;
enddo ;
enddo
305 evap_cfl_limit, minimum_forcing_depth)
321 type(
surface),
intent(inout) :: state
323 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
331 integer :: m, is, ie, js, je, isd, ied, jsd, jed
332 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
333 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
335 if (.not.
associated(cs))
return
337 if (cs%coupled_tracers)
then
341 call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
342 state%tr_fields, idim=(/isd, is, ie, ied/), &
343 jdim=(/jsd, js, je, jed/) )
353 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
354 real,
dimension(:),
intent(out) :: stocks
359 character(len=*),
dimension(:),
intent(out) :: names
360 character(len=*),
dimension(:),
intent(out) :: units
361 integer,
optional,
intent(in) :: stock_index
364 integer :: i, j, k, is, ie, js, je, nz, m
365 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
368 if (.not.
associated(cs))
return
369 if (cs%ntr < 1)
return
371 if (
present(stock_index))
then ;
if (stock_index > 0)
then
379 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller=
"advection_test_stock")
381 do k=1,nz ;
do j=js,je ;
do i=is,ie
382 stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
383 (g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j) * h(i,j,k))
384 enddo ;
enddo ;
enddo
385 stocks(m) = gv%H_to_kg_m2 * stocks(m)
397 if (
associated(cs))
then
398 if (
associated(cs%tr))
deallocate(cs%tr)