24 use coupler_types_mod,
only : coupler_type_set_data, ind_csurf
27 implicit none ;
private
29 #include <MOM_memory.h>
31 public register_oil_tracer, initialize_oil_tracer
32 public oil_tracer_column_physics, oil_tracer_surface_state
33 public oil_stock, oil_tracer_end
35 integer,
parameter :: ntr_max = 20
40 logical :: coupled_tracers = .false.
41 character(len=200) :: ic_file
44 real :: oil_source_longitude
45 real :: oil_source_latitude
46 integer :: oil_source_i=-999
47 integer :: oil_source_j=-999
48 real :: oil_source_rate
49 real :: oil_start_year
53 type(time_type),
pointer :: time => null()
55 real,
pointer :: tr(:,:,:,:) => null()
56 real,
dimension(NTR_MAX) :: ic_val = 0.0
57 real,
dimension(NTR_MAX) :: young_val = 0.0
58 real,
dimension(NTR_MAX) :: land_val = -1.0
59 real,
dimension(NTR_MAX) :: sfc_growth_rate
60 real,
dimension(NTR_MAX) :: oil_decay_days
61 real,
dimension(NTR_MAX) :: oil_decay_rate
62 integer,
dimension(NTR_MAX) :: oil_source_k
63 logical :: oil_may_reinit
65 integer,
dimension(NTR_MAX) :: ind_tr
77 function register_oil_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
89 character(len=40) :: mdl =
"oil_tracer"
91 #include "version_variable.h"
92 character(len=200) :: inputdir
93 character(len=48) :: var_name
94 character(len=3) :: name_tag
95 character(len=48) :: flux_units
97 real,
pointer :: tr_ptr(:,:,:) => null()
98 logical :: register_oil_tracer
99 integer :: isd, ied, jsd, jed, nz, m, i, j
100 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
102 if (
associated(cs))
then
103 call mom_error(warning,
"register_oil_tracer called with an "// &
104 "associated control structure.")
111 call get_param(param_file, mdl,
"OIL_IC_FILE", cs%IC_file, &
112 "The file in which the oil tracer initial values can be "//&
113 "found, or an empty string for internal initialization.", &
115 if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file,
'/') == 0))
then
117 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
118 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
119 call log_param(param_file, mdl,
"INPUTDIR/CFC_IC_FILE", cs%IC_file)
121 call get_param(param_file, mdl,
"OIL_IC_FILE_IS_Z", cs%Z_IC_file, &
122 "If true, OIL_IC_FILE is in depth space, not layer space", &
125 call get_param(param_file, mdl,
"OIL_MAY_REINIT", cs%oil_may_reinit, &
126 "If true, oil tracers may go through the initialization "//&
127 "code if they are not found in the restart files. "//&
128 "Otherwise it is a fatal error if the oil tracers are not "//&
129 "found in the restart files of a restarted run.", &
131 call get_param(param_file, mdl,
"OIL_SOURCE_LONGITUDE", cs%oil_source_longitude, &
132 "The geographic longitude of the oil source.", units=
"degrees E", &
133 fail_if_missing=.true.)
134 call get_param(param_file, mdl,
"OIL_SOURCE_LATITUDE", cs%oil_source_latitude, &
135 "The geographic latitude of the oil source.", units=
"degrees N", &
136 fail_if_missing=.true.)
137 call get_param(param_file, mdl,
"OIL_SOURCE_LAYER", cs%oil_source_k, &
138 "The layer into which the oil is introduced, or a "//&
139 "negative number for a vertically uniform source, "//&
140 "or 0 not to use this tracer.", units=
"Layer", default=0)
141 call get_param(param_file, mdl,
"OIL_SOURCE_RATE", cs%oil_source_rate, &
142 "The rate of oil injection.", units=
"kg s-1", default=1.0)
143 call get_param(param_file, mdl,
"OIL_DECAY_DAYS", cs%oil_decay_days, &
144 "The decay timescale in days (if positive), or no decay "//&
145 "if 0, or use the temperature dependent decay rate of "//&
146 "Adcroft et al. (GRL, 2010) if negative.", units=
"days", &
148 call get_param(param_file, mdl,
"OIL_DATED_START_YEAR", cs%oil_start_year, &
149 "The time at which the oil source starts", units=
"years", &
151 call get_param(param_file, mdl,
"OIL_DATED_END_YEAR", cs%oil_end_year, &
152 "The time at which the oil source ends", units=
"years", &
156 cs%oil_decay_rate(:) = 0.
158 if (cs%oil_source_k(m)/=0)
then
159 write(name_tag(1:3),
'("_",I2.2)') m
161 cs%tr_desc(m) = var_desc(
"oil"//trim(name_tag),
"kg m-3",
"Oil Tracer", caller=mdl)
163 if (cs%oil_decay_days(m)>0.)
then
164 cs%oil_decay_rate(m)=1./(86400.0*cs%oil_decay_days(m))
165 elseif (cs%oil_decay_days(m)<0.)
then
166 cs%oil_decay_rate(m)=-1.
170 call log_param(param_file, mdl,
"OIL_DECAY_RATE", cs%oil_decay_rate(1:cs%ntr))
173 if (gv%Boussinesq)
then ; flux_units =
"kg s-1"
174 else ; flux_units =
"kg m-3 kg s-1" ;
endif
176 allocate(cs%tr(isd:ied,jsd:jed,nz,cs%ntr)) ; cs%tr(:,:,:,:) = 0.0
181 tr_ptr => cs%tr(:,:,:,m)
182 call query_vardesc(cs%tr_desc(m), name=var_name, caller=
"register_oil_tracer")
184 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, tr_desc=cs%tr_desc(m), &
185 registry_diags=.true., flux_units=flux_units, restart_cs=restart_cs, &
186 mandatory=.not.cs%oil_may_reinit)
191 if (cs%coupled_tracers) &
192 cs%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//
'_flux', &
193 flux_type=
' ', implementation=
' ', caller=
"register_oil_tracer")
197 cs%restart_CSp => restart_cs
198 register_oil_tracer = .true.
200 end function register_oil_tracer
203 subroutine initialize_oil_tracer(restart, day, G, GV, US, h, diag, OBC, CS, &
205 logical,
intent(in) :: restart
207 type(time_type),
target,
intent(in) :: day
211 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
213 type(
diag_ctrl),
target,
intent(in) :: diag
223 character(len=16) :: name
224 character(len=72) :: longname
225 character(len=48) :: units
226 character(len=48) :: flux_units
229 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
230 integer :: isdb, iedb, jsdb, jedb
232 if (.not.
associated(cs))
return
233 if (cs%ntr < 1)
return
234 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
235 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
236 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
239 do j=g%jsdB+1,g%jed ;
do i=g%isdB+1,g%ied
242 if (cs%oil_source_longitude<g%geoLonBu(i,j) .and. &
243 cs%oil_source_longitude>=g%geoLonBu(i-1,j) .and. &
244 cs%oil_source_latitude<g%geoLatBu(i,j) .and. &
245 cs%oil_source_latitude>=g%geoLatBu(i,j-1) )
then
255 call query_vardesc(cs%tr_desc(m), name=name, caller=
"initialize_oil_tracer")
256 if ((.not.restart) .or. (cs%oil_may_reinit .and. .not. &
259 if (len_trim(cs%IC_file) > 0)
then
262 call mom_error(fatal,
"initialize_oil_tracer: "// &
263 "Unable to open "//cs%IC_file)
265 if (cs%Z_IC_file)
then
266 ok = tracer_z_init(cs%tr(:,:,:,m), h, cs%IC_file, name, &
269 ok = tracer_z_init(cs%tr(:,:,:,m), h, cs%IC_file, &
270 trim(name), g, us, -1e34, 0.0)
271 if (.not.ok)
call mom_error(fatal,
"initialize_oil_tracer: "//&
272 "Unable to read "//trim(name)//
" from "//&
273 trim(cs%IC_file)//
".")
276 call mom_read_data(cs%IC_file, trim(name), cs%tr(:,:,:,m), g%Domain)
279 do k=1,nz ;
do j=js,je ;
do i=is,ie
280 if (g%mask2dT(i,j) < 0.5)
then
281 cs%tr(i,j,k,m) = cs%land_val(m)
283 cs%tr(i,j,k,m) = cs%IC_val(m)
285 enddo ;
enddo ;
enddo
291 if (
associated(obc))
then
295 end subroutine initialize_oil_tracer
298 subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, &
299 evap_CFL_limit, minimum_forcing_depth)
302 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
304 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
306 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
310 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
314 type(
forcing),
intent(in) :: fluxes
316 real,
intent(in) :: dt
320 real,
optional,
intent(in) :: evap_cfl_limit
322 real,
optional,
intent(in) :: minimum_forcing_depth
332 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
333 real :: isecs_per_year = 1.0 / (365.0*86400.0)
334 real :: year, h_total, ldecay
335 integer :: i, j, k, is, ie, js, je, nz, m, k_max
336 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
338 if (.not.
associated(cs))
return
339 if (cs%ntr < 1)
return
341 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
343 do k=1,nz ;
do j=js,je ;
do i=is,ie
344 h_work(i,j,k) = h_old(i,j,k)
345 enddo ;
enddo ;
enddo
346 call applytracerboundaryfluxesinout(g, gv, cs%tr(:,:,:,m) , dt, fluxes, h_work, &
347 evap_cfl_limit, minimum_forcing_depth)
348 call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
352 call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
356 year = time_type_to_real(cs%Time) * isecs_per_year
360 do k=1,nz ;
do j=js,je ;
do i=is,ie
363 if (cs%oil_decay_rate(m)>0.)
then
364 cs%tr(i,j,k,m) = g%mask2dT(i,j)*max(1.-dt*cs%oil_decay_rate(m),0.)*cs%tr(i,j,k,m)
365 elseif (cs%oil_decay_rate(m)<0.)
then
366 ldecay = 12.*(3.0**(-(tv%T(i,j,k)-20.)/10.))
367 ldecay = 1./(86400.*ldecay)
368 cs%tr(i,j,k,m) = g%mask2dT(i,j)*max(1.-dt*ldecay,0.)*cs%tr(i,j,k,m)
370 enddo ;
enddo ;
enddo
374 if (year>=cs%oil_start_year .and. year<=cs%oil_end_year .and. &
375 cs%oil_source_i>-999 .and. cs%oil_source_j>-999)
then
376 i=cs%oil_source_i ; j=cs%oil_source_j
377 k_max=nz ; h_total=0.
379 h_total = h_total + h_new(i,j,k)
380 if (h_total<10.) k_max=k-1
386 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + cs%oil_source_rate*dt / &
387 ((h_new(i,j,k)+gv%H_subroundoff) * g%areaT(i,j) )
389 h_total=gv%H_subroundoff
391 h_total = h_total + h_new(i,j,k)
394 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + cs%oil_source_rate*dt/(h_total &
401 end subroutine oil_tracer_column_physics
405 function oil_stock(h, stocks, G, GV, CS, names, units, stock_index)
408 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
409 real,
dimension(:),
intent(out) :: stocks
413 character(len=*),
dimension(:),
intent(out) :: names
414 character(len=*),
dimension(:),
intent(out) :: units
415 integer,
optional,
intent(in) :: stock_index
424 integer :: i, j, k, is, ie, js, je, nz, m
425 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
428 if (.not.
associated(cs))
return
429 if (cs%ntr < 1)
return
431 if (
present(stock_index))
then ;
if (stock_index > 0)
then
439 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller=
"oil_stock")
440 units(m) = trim(units(m))//
" kg"
442 do k=1,nz ;
do j=js,je ;
do i=is,ie
443 stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
444 (g%mask2dT(i,j) * g%areaT(i,j) * h(i,j,k))
445 enddo ;
enddo ;
enddo
446 stocks(m) = gv%H_to_kg_m2 * stocks(m)
450 end function oil_stock
455 subroutine oil_tracer_surface_state(state, h, G, CS)
457 type(
surface),
intent(inout) :: state
459 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
467 integer :: m, is, ie, js, je, isd, ied, jsd, jed
468 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
469 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
471 if (.not.
associated(cs))
return
473 if (cs%coupled_tracers)
then
477 call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
478 state%tr_fields, idim=(/isd, is, ie, ied/), &
479 jdim=(/jsd, js, je, jed/) )
483 end subroutine oil_tracer_surface_state
486 subroutine oil_tracer_end(CS)
491 if (
associated(cs))
then
492 if (
associated(cs%tr))
deallocate(cs%tr)
495 end subroutine oil_tracer_end