24 use coupler_types_mod,
only : coupler_type_set_data, ind_csurf
27 implicit none ;
private
29 #include <MOM_memory.h>
31 public register_boundary_impulse_tracer, initialize_boundary_impulse_tracer
32 public boundary_impulse_tracer_column_physics, boundary_impulse_tracer_surface_state
33 public boundary_impulse_stock, boundary_impulse_tracer_end
36 integer,
parameter :: ntr_max = 1
40 integer :: ntr=ntr_max
41 logical :: coupled_tracers = .false.
42 type(time_type),
pointer :: time => null()
44 real,
pointer :: tr(:,:,:,:) => null()
45 logical :: tracers_may_reinit
46 integer,
dimension(NTR_MAX) :: ind_tr
50 real,
dimension(NTR_MAX) :: land_val = -1.0
52 real :: remaining_source_time
65 function register_boundary_impulse_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
77 character(len=40) :: mdl =
"boundary_impulse_tracer"
78 character(len=200) :: inputdir
79 character(len=48) :: var_name
80 character(len=3) :: name_tag
81 character(len=48) :: flux_units
84 #include "version_variable.h"
85 real,
pointer :: tr_ptr(:,:,:) => null()
86 real,
pointer :: rem_time_ptr => null()
87 logical :: register_boundary_impulse_tracer
88 integer :: isd, ied, jsd, jed, nz, m, i, j
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_boundary_impulse_tracer called with an "// &
93 "associated control structure.")
100 call get_param(param_file, mdl,
"IMPULSE_SOURCE_TIME", cs%remaining_source_time, &
101 "Length of time for the boundary tracer to be injected "//&
102 "into the mixed layer. After this time has elapsed, the "//&
103 "surface becomes a sink for the boundary impulse tracer.", &
105 call get_param(param_file, mdl,
"TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
106 "If true, tracers may go through the initialization code "//&
107 "if they are not found in the restart files. Otherwise "//&
108 "it is a fatal error if the tracers are not found in the "//&
109 "restart files of a restarted run.", default=.false.)
111 allocate(cs%tr(isd:ied,jsd:jed,nz,cs%ntr)) ; cs%tr(:,:,:,:) = 0.0
113 cs%nkml = max(gv%nkml,1)
118 cs%tr_desc(m) = var_desc(trim(
"boundary_impulse"),
"kg kg-1", &
119 "Boundary impulse tracer", caller=mdl)
120 if (gv%Boussinesq)
then ; flux_units =
"kg kg-1 m3 s-1"
121 else ; flux_units =
"kg s-1" ;
endif
123 tr_ptr => cs%tr(:,:,:,m)
124 call query_vardesc(cs%tr_desc(m), name=var_name, caller=
"register_boundary_impulse_tracer")
126 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, tr_desc=cs%tr_desc(m), &
127 registry_diags=.true., flux_units=flux_units, &
128 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
133 if (cs%coupled_tracers) &
134 cs%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//
'_flux', &
135 flux_type=
' ', implementation=
' ', caller=
"register_boundary_impulse_tracer")
138 rem_time_ptr => cs%remaining_source_time
140 .not.cs%tracers_may_reinit, restart_cs, &
141 "Remaining time to apply BIR source",
"s")
144 cs%restart_CSp => restart_cs
145 register_boundary_impulse_tracer = .true.
147 end function register_boundary_impulse_tracer
150 subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, h, diag, OBC, CS, &
152 logical,
intent(in) :: restart
154 type(time_type),
target,
intent(in) :: day
157 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
159 type(
diag_ctrl),
target,
intent(in) :: diag
170 character(len=16) :: name
171 character(len=72) :: longname
172 character(len=48) :: units
173 character(len=48) :: flux_units
176 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
177 integer :: isdb, iedb, jsdb, jedb
179 if (.not.
associated(cs))
return
180 if (cs%ntr < 1)
return
181 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
182 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
183 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
187 name =
"boundary_impulse"
190 call query_vardesc(cs%tr_desc(m), name=name, caller=
"initialize_boundary_impulse_tracer")
191 if ((.not.restart) .or. (.not. &
193 do k=1,cs%nkml ;
do j=jsd,jed ;
do i=isd,ied
195 enddo ;
enddo ;
enddo
199 if (
associated(obc))
then
203 end subroutine initialize_boundary_impulse_tracer
206 subroutine boundary_impulse_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
207 tv, debug, evap_CFL_limit, minimum_forcing_depth)
210 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
212 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
214 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
218 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
222 type(
forcing),
intent(in) :: fluxes
224 real,
intent(in) :: dt
229 logical,
intent(in) :: debug
230 real,
optional,
intent(in) :: evap_cfl_limit
232 real,
optional,
intent(in) :: minimum_forcing_depth
243 real :: isecs_per_year = 1.0 / (365.0*86400.0)
244 real :: year, h_total, scale, htot, ih_limit
245 integer :: secs, days
246 integer :: i, j, k, is, ie, js, je, nz, m, k_max
247 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
249 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
251 if (.not.
associated(cs))
return
252 if (cs%ntr < 1)
return
255 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
256 do k=1,nz ;
do j=js,je ;
do i=is,ie
257 h_work(i,j,k) = h_old(i,j,k)
258 enddo ;
enddo ;
enddo
259 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,1), dt, fluxes, h_work, &
260 evap_cfl_limit, minimum_forcing_depth)
261 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,1), g, gv)
263 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,1), g, gv)
268 if (cs%remaining_source_time>0.0)
then
269 do k=1,cs%nkml ;
do j=js,je ;
do i=is,ie
271 enddo ;
enddo ;
enddo
272 cs%remaining_source_time = cs%remaining_source_time-dt
274 do k=1,cs%nkml ;
do j=js,je ;
do i=is,ie
276 enddo ;
enddo ;
enddo
281 end subroutine boundary_impulse_tracer_column_physics
284 function boundary_impulse_stock(h, stocks, G, GV, CS, names, units, stock_index)
287 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in ) :: h
288 real,
dimension(:),
intent( out) :: stocks
292 character(len=*),
dimension(:),
intent( out) :: names
293 character(len=*),
dimension(:),
intent( out) :: units
294 integer,
optional,
intent(in ) :: stock_index
296 integer :: boundary_impulse_stock
303 integer :: i, j, k, is, ie, js, je, nz, m
304 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
306 boundary_impulse_stock = 0
307 if (.not.
associated(cs))
return
308 if (cs%ntr < 1)
return
310 if (
present(stock_index))
then ;
if (stock_index > 0)
then
318 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller=
"boundary_impulse_stock")
319 units(m) = trim(units(m))//
" kg"
321 do k=1,nz ;
do j=js,je ;
do i=is,ie
322 stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
323 (g%mask2dT(i,j) * g%areaT(i,j) * h(i,j,k))
324 enddo ;
enddo ;
enddo
325 stocks(m) = gv%H_to_kg_m2 * stocks(m)
328 boundary_impulse_stock = cs%ntr
330 end function boundary_impulse_stock
335 subroutine boundary_impulse_tracer_surface_state(state, h, G, CS)
337 type(
surface),
intent(inout) :: state
339 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
347 integer :: m, is, ie, js, je, isd, ied, jsd, jed
348 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
349 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
351 if (.not.
associated(cs))
return
353 if (cs%coupled_tracers)
then
357 call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
358 state%tr_fields, idim=(/isd, is, ie, ied/), &
359 jdim=(/jsd, js, je, jed/) )
363 end subroutine boundary_impulse_tracer_surface_state
366 subroutine boundary_impulse_tracer_end(CS)
371 if (
associated(cs))
then
372 if (
associated(cs%tr))
deallocate(cs%tr)
375 end subroutine boundary_impulse_tracer_end