22 use coupler_types_mod,
only : coupler_type_set_data, ind_csurf
25 implicit none ;
private
27 #include <MOM_memory.h>
29 public register_advection_test_tracer, initialize_advection_test_tracer
30 public advection_test_tracer_surface_state, advection_test_tracer_end
31 public advection_test_tracer_column_physics, advection_test_stock
33 integer,
parameter :: ntr = 11
38 logical :: coupled_tracers = .false.
39 character(len=200) :: tracer_ic_file
40 type(time_type),
pointer :: time => null()
42 real,
pointer :: tr(:,:,:,:) => null()
43 real :: land_val(ntr) = -1.0
45 logical :: tracers_may_reinit
53 integer,
dimension(NTR) :: ind_tr
66 function register_advection_test_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
78 character(len=80) :: name, longname
80 #include "version_variable.h"
81 character(len=40) :: mdl =
"advection_test_tracer"
82 character(len=200) :: inputdir
83 character(len=48) :: flux_units
85 real,
pointer :: tr_ptr(:,:,:) => null()
86 logical :: register_advection_test_tracer
87 integer :: isd, ied, jsd, jed, nz, m
88 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
90 if (
associated(cs))
then
91 call mom_error(warning,
"register_advection_test_tracer called with an "// &
92 "associated control structure.")
100 call get_param(param_file, mdl,
"ADVECTION_TEST_X_ORIGIN", cs%x_origin, &
101 "The x-coorindate of the center of the test-functions.", default=0.)
102 call get_param(param_file, mdl,
"ADVECTION_TEST_Y_ORIGIN", cs%y_origin, &
103 "The y-coorindate of the center of the test-functions.", default=0.)
104 call get_param(param_file, mdl,
"ADVECTION_TEST_X_WIDTH", cs%x_width, &
105 "The x-width of the test-functions.", default=0.)
106 call get_param(param_file, mdl,
"ADVECTION_TEST_Y_WIDTH", cs%y_width, &
107 "The y-width of the test-functions.", default=0.)
108 call get_param(param_file, mdl,
"ADVECTION_TEST_TRACER_IC_FILE", cs%tracer_IC_file, &
109 "The name of a file from which to read the initial "//&
110 "conditions for the tracers, or blank to initialize "//&
111 "them internally.", default=
" ")
113 if (len_trim(cs%tracer_IC_file) >= 1)
then
114 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
115 cs%tracer_IC_file = trim(slasher(inputdir))//trim(cs%tracer_IC_file)
116 call log_param(param_file, mdl,
"INPUTDIR/ADVECTION_TEST_TRACER_IC_FILE", &
119 call get_param(param_file, mdl,
"SPONGE", cs%use_sponge, &
120 "If true, sponges may be applied anywhere in the domain. "//&
121 "The exact location and properties of those sponges are "//&
122 "specified from MOM_initialization.F90.", default=.false.)
124 call get_param(param_file, mdl,
"TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
125 "If true, tracers may go through the initialization code "//&
126 "if they are not found in the restart files. Otherwise "//&
127 "it is a fatal error if the tracers are not found in the "//&
128 "restart files of a restarted run.", default=.false.)
131 allocate(cs%tr(isd:ied,jsd:jed,nz,ntr)) ; cs%tr(:,:,:,:) = 0.0
134 if (m < 10)
then ;
write(name,
'("tr",I1.1)') m
135 else ;
write(name,
'("tr",I2.2)') m ;
endif
136 write(longname,
'("Concentration of Tracer ",I2.2)') m
137 cs%tr_desc(m) = var_desc(name, units=
"kg kg-1", longname=longname, caller=mdl)
138 if (gv%Boussinesq)
then ; flux_units =
"kg kg-1 m3 s-1"
139 else ; flux_units =
"kg s-1" ;
endif
144 tr_ptr => cs%tr(:,:,:,m)
146 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
147 name=name, longname=longname, units=
"kg kg-1", &
148 registry_diags=.true., flux_units=flux_units, &
149 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
154 if (cs%coupled_tracers) &
155 cs%ind_tr(m) = aof_set_coupler_flux(trim(name)//
'_flux', &
156 flux_type=
' ', implementation=
' ', caller=
"register_advection_test_tracer")
160 cs%restart_CSp => restart_cs
161 register_advection_test_tracer = .true.
162 end function register_advection_test_tracer
165 subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS, &
167 logical,
intent(in) :: restart
169 type(time_type),
target,
intent(in) :: day
172 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
174 type(
diag_ctrl),
target,
intent(in) :: diag
184 real,
allocatable :: temp(:,:,:)
185 real,
pointer,
dimension(:,:,:) :: &
186 obc_tr1_u => null(), &
190 character(len=16) :: name
191 character(len=72) :: longname
192 character(len=48) :: units
193 character(len=48) :: flux_units
195 real,
pointer :: tr_ptr(:,:,:) => null()
201 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
202 integer :: isdb, iedb, jsdb, jedb
203 real :: tmpx, tmpy, locx, locy
205 if (.not.
associated(cs))
return
206 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
207 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
208 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
209 h_neglect = gv%H_subroundoff
214 call query_vardesc(cs%tr_desc(m), name=name, &
215 caller=
"initialize_advection_test_tracer")
216 if ((.not.restart) .or. (cs%tracers_may_reinit .and. .not. &
218 do k=1,nz ;
do j=js,je ;
do i=is,ie
220 enddo ;
enddo ;
enddo
222 do j=js,je ;
do i=is,ie
223 if (abs(g%geoLonT(i,j)-cs%x_origin)<0.5*cs%x_width .and. &
224 abs(g%geoLatT(i,j)-cs%y_origin)<0.5*cs%y_width) cs%tr(i,j,k,m) = 1.0
227 do j=js,je ;
do i=is,ie
228 locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
229 locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
230 cs%tr(i,j,k,m) = max(0.0, 1.0-locx)*max(0.0, 1.0-locy)
233 do j=js,je ;
do i=is,ie
234 locx=min(1.0, abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width)*(acos(0.0)*2.)
235 locy=min(1.0, abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width)*(acos(0.0)*2.)
236 cs%tr(i,j,k,m) = (1.0+cos(locx))*(1.0+cos(locy))*0.25
239 do j=js,je ;
do i=is,ie
240 locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
241 locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
242 if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
245 do j=js,je ;
do i=is,ie
246 locx=(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
247 locy=(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
248 if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
249 if (locx>0.0.and.abs(locy)<0.2) cs%tr(i,j,k,m) = 0.0
255 end subroutine initialize_advection_test_tracer
260 subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
261 evap_CFL_limit, minimum_forcing_depth)
264 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
266 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
268 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
272 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
276 type(
forcing),
intent(in) :: fluxes
278 real,
intent(in) :: dt
281 real,
optional,
intent(in) :: evap_cfl_limit
283 real,
optional,
intent(in) :: minimum_forcing_depth
292 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
294 real :: c1(szi_(g),szk_(g))
295 integer :: i, j, k, is, ie, js, je, nz, m
296 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
298 if (.not.
associated(cs))
return
300 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
302 do k=1,nz ;
do j=js,je ;
do i=is,ie
303 h_work(i,j,k) = h_old(i,j,k)
304 enddo ;
enddo ;
enddo
305 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m) , dt, fluxes, h_work, &
306 evap_cfl_limit, minimum_forcing_depth)
307 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
311 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
315 end subroutine advection_test_tracer_column_physics
320 subroutine advection_test_tracer_surface_state(state, h, G, CS)
322 type(
surface),
intent(inout) :: state
324 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
332 integer :: m, is, ie, js, je, isd, ied, jsd, jed
333 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
334 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
336 if (.not.
associated(cs))
return
338 if (cs%coupled_tracers)
then
342 call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
343 state%tr_fields, idim=(/isd, is, ie, ied/), &
344 jdim=(/jsd, js, je, jed/) )
348 end subroutine advection_test_tracer_surface_state
352 function advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
354 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
355 real,
dimension(:),
intent(out) :: stocks
360 character(len=*),
dimension(:),
intent(out) :: names
361 character(len=*),
dimension(:),
intent(out) :: units
362 integer,
optional,
intent(in) :: stock_index
363 integer :: advection_test_stock
365 integer :: i, j, k, is, ie, js, je, nz, m
366 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
368 advection_test_stock = 0
369 if (.not.
associated(cs))
return
370 if (cs%ntr < 1)
return
372 if (
present(stock_index))
then ;
if (stock_index > 0)
then
380 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller=
"advection_test_stock")
382 do k=1,nz ;
do j=js,je ;
do i=is,ie
383 stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
384 (g%mask2dT(i,j) * g%areaT(i,j) * h(i,j,k))
385 enddo ;
enddo ;
enddo
386 stocks(m) = gv%H_to_kg_m2 * stocks(m)
388 advection_test_stock = cs%ntr
390 end function advection_test_stock
393 subroutine advection_test_tracer_end(CS)
398 if (
associated(cs))
then
399 if (
associated(cs%tr))
deallocate(cs%tr)
402 end subroutine advection_test_tracer_end