24 use coupler_types_mod,
only : coupler_type_set_data, ind_csurf
27 implicit none ;
private
29 #include <MOM_memory.h>
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
90 character(len=40) :: mdl =
"oil_tracer"
92 #include "version_variable.h"
93 character(len=200) :: inputdir
94 character(len=48) :: var_name
95 character(len=3) :: name_tag
96 character(len=48) :: flux_units
98 real,
pointer :: tr_ptr(:,:,:) => null()
100 integer :: isd, ied, jsd, jed, nz, m, i, j
101 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
103 if (
associated(cs))
then
104 call mom_error(warning,
"register_oil_tracer called with an "// &
105 "associated control structure.")
112 call get_param(param_file, mdl,
"OIL_IC_FILE", cs%IC_file, &
113 "The file in which the oil tracer initial values can be "//&
114 "found, or an empty string for internal initialization.", &
116 if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file,
'/') == 0))
then
118 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
119 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
120 call log_param(param_file, mdl,
"INPUTDIR/CFC_IC_FILE", cs%IC_file)
122 call get_param(param_file, mdl,
"OIL_IC_FILE_IS_Z", cs%Z_IC_file, &
123 "If true, OIL_IC_FILE is in depth space, not layer space", &
126 call get_param(param_file, mdl,
"OIL_MAY_REINIT", cs%oil_may_reinit, &
127 "If true, oil tracers may go through the initialization "//&
128 "code if they are not found in the restart files. "//&
129 "Otherwise it is a fatal error if the oil tracers are not "//&
130 "found in the restart files of a restarted run.", &
132 call get_param(param_file, mdl,
"OIL_SOURCE_LONGITUDE", cs%oil_source_longitude, &
133 "The geographic longitude of the oil source.", units=
"degrees E", &
134 fail_if_missing=.true.)
135 call get_param(param_file, mdl,
"OIL_SOURCE_LATITUDE", cs%oil_source_latitude, &
136 "The geographic latitude of the oil source.", units=
"degrees N", &
137 fail_if_missing=.true.)
138 call get_param(param_file, mdl,
"OIL_SOURCE_LAYER", cs%oil_source_k, &
139 "The layer into which the oil is introduced, or a "//&
140 "negative number for a vertically uniform source, "//&
141 "or 0 not to use this tracer.", units=
"Layer", default=0)
142 call get_param(param_file, mdl,
"OIL_SOURCE_RATE", cs%oil_source_rate, &
143 "The rate of oil injection.", units=
"kg s-1", scale=us%T_to_s, default=1.0)
144 call get_param(param_file, mdl,
"OIL_DECAY_DAYS", cs%oil_decay_days, &
145 "The decay timescale in days (if positive), or no decay "//&
146 "if 0, or use the temperature dependent decay rate of "//&
147 "Adcroft et al. (GRL, 2010) if negative.", units=
"days", &
149 call get_param(param_file, mdl,
"OIL_DATED_START_YEAR", cs%oil_start_year, &
150 "The time at which the oil source starts", units=
"years", &
152 call get_param(param_file, mdl,
"OIL_DATED_END_YEAR", cs%oil_end_year, &
153 "The time at which the oil source ends", units=
"years", &
157 cs%oil_decay_rate(:) = 0.
159 if (cs%oil_source_k(m)/=0)
then
160 write(name_tag(1:3),
'("_",I2.2)') m
162 cs%tr_desc(m) =
var_desc(
"oil"//trim(name_tag),
"kg m-3",
"Oil Tracer", caller=mdl)
164 if (cs%oil_decay_days(m)>0.)
then
165 cs%oil_decay_rate(m) = 1. / (86400.0*us%s_to_T * cs%oil_decay_days(m))
166 elseif (cs%oil_decay_days(m)<0.)
then
167 cs%oil_decay_rate(m) = -1.
171 call log_param(param_file, mdl,
"OIL_DECAY_RATE", us%s_to_T*cs%oil_decay_rate(1:cs%ntr))
174 if (gv%Boussinesq)
then ; flux_units =
"kg s-1"
175 else ; flux_units =
"kg m-3 kg s-1" ;
endif
177 allocate(cs%tr(isd:ied,jsd:jed,nz,cs%ntr)) ; cs%tr(:,:,:,:) = 0.0
182 tr_ptr => cs%tr(:,:,:,m)
183 call query_vardesc(cs%tr_desc(m), name=var_name, caller=
"register_oil_tracer")
185 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, tr_desc=cs%tr_desc(m), &
186 registry_diags=.true., flux_units=flux_units, restart_cs=restart_cs, &
187 mandatory=.not.cs%oil_may_reinit)
192 if (cs%coupled_tracers) &
194 flux_type=
' ', implementation=
' ', caller=
"register_oil_tracer")
198 cs%restart_CSp => restart_cs
206 logical,
intent(in) :: restart
208 type(time_type),
target,
intent(in) :: day
212 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
214 type(
diag_ctrl),
target,
intent(in) :: diag
224 character(len=16) :: name
225 character(len=72) :: longname
226 character(len=48) :: units
227 character(len=48) :: flux_units
230 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
231 integer :: isdb, iedb, jsdb, jedb
233 if (.not.
associated(cs))
return
234 if (cs%ntr < 1)
return
235 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
236 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
237 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
240 do j=g%jsdB+1,g%jed ;
do i=g%isdB+1,g%ied
243 if (cs%oil_source_longitude<g%geoLonBu(i,j) .and. &
244 cs%oil_source_longitude>=g%geoLonBu(i-1,j) .and. &
245 cs%oil_source_latitude<g%geoLatBu(i,j) .and. &
246 cs%oil_source_latitude>=g%geoLatBu(i,j-1) )
then
256 call query_vardesc(cs%tr_desc(m), name=name, caller=
"initialize_oil_tracer")
257 if ((.not.restart) .or. (cs%oil_may_reinit .and. .not. &
260 if (len_trim(cs%IC_file) > 0)
then
263 call mom_error(fatal,
"initialize_oil_tracer: "// &
264 "Unable to open "//cs%IC_file)
266 if (cs%Z_IC_file)
then
271 trim(name), g, us, -1e34, 0.0)
272 if (.not.ok)
call mom_error(fatal,
"initialize_oil_tracer: "//&
273 "Unable to read "//trim(name)//
" from "//&
274 trim(cs%IC_file)//
".")
277 call mom_read_data(cs%IC_file, trim(name), cs%tr(:,:,:,m), g%Domain)
280 do k=1,nz ;
do j=js,je ;
do i=is,ie
281 if (g%mask2dT(i,j) < 0.5)
then
282 cs%tr(i,j,k,m) = cs%land_val(m)
284 cs%tr(i,j,k,m) = cs%IC_val(m)
286 enddo ;
enddo ;
enddo
292 if (
associated(obc))
then
299 subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, tv, &
300 evap_CFL_limit, minimum_forcing_depth)
303 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
305 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
307 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
311 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
315 type(
forcing),
intent(in) :: fluxes
317 real,
intent(in) :: dt
322 real,
optional,
intent(in) :: evap_cfl_limit
324 real,
optional,
intent(in) :: minimum_forcing_depth
334 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
335 real :: isecs_per_year = 1.0 / (365.0*86400.0)
336 real :: year, h_total, ldecay
337 integer :: i, j, k, is, ie, js, je, nz, m, k_max
338 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
340 if (.not.
associated(cs))
return
341 if (cs%ntr < 1)
return
343 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
345 do k=1,nz ;
do j=js,je ;
do i=is,ie
346 h_work(i,j,k) = h_old(i,j,k)
347 enddo ;
enddo ;
enddo
349 evap_cfl_limit, minimum_forcing_depth)
358 year = time_type_to_real(cs%Time) * isecs_per_year
362 do k=1,nz ;
do j=js,je ;
do i=is,ie
365 if (cs%oil_decay_rate(m)>0.)
then
366 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)
367 elseif (cs%oil_decay_rate(m)<0.)
then
368 ldecay = 12.*(3.0**(-(tv%T(i,j,k)-20.)/10.))
369 ldecay = 1. / (86400.*us%s_to_T * ldecay)
370 cs%tr(i,j,k,m) = g%mask2dT(i,j)*max(1. - dt*ldecay,0.)*cs%tr(i,j,k,m)
372 enddo ;
enddo ;
enddo
376 if (year>=cs%oil_start_year .and. year<=cs%oil_end_year .and. &
377 cs%oil_source_i>-999 .and. cs%oil_source_j>-999)
then
378 i=cs%oil_source_i ; j=cs%oil_source_j
379 k_max=nz ; h_total=0.
381 h_total = h_total + h_new(i,j,k)
382 if (h_total<10.) k_max=k-1
388 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + cs%oil_source_rate*dt / &
389 ((h_new(i,j,k)+gv%H_subroundoff) * g%US%L_to_m**2*g%areaT(i,j) )
391 h_total=gv%H_subroundoff
393 h_total = h_total + h_new(i,j,k)
396 cs%tr(i,j,k,m) = cs%tr(i,j,k,m) + cs%oil_source_rate*dt/(h_total &
397 * g%US%L_to_m**2*g%areaT(i,j) )
407 function oil_stock(h, stocks, G, GV, CS, names, units, stock_index)
410 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
411 real,
dimension(:),
intent(out) :: stocks
415 character(len=*),
dimension(:),
intent(out) :: names
416 character(len=*),
dimension(:),
intent(out) :: units
417 integer,
optional,
intent(in) :: stock_index
426 integer :: i, j, k, is, ie, js, je, nz, m
427 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
430 if (.not.
associated(cs))
return
431 if (cs%ntr < 1)
return
433 if (
present(stock_index))
then ;
if (stock_index > 0)
then
441 call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller=
"oil_stock")
442 units(m) = trim(units(m))//
" kg"
444 do k=1,nz ;
do j=js,je ;
do i=is,ie
445 stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
446 (g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j) * h(i,j,k))
447 enddo ;
enddo ;
enddo
448 stocks(m) = gv%H_to_kg_m2 * stocks(m)
459 type(
surface),
intent(inout) :: state
461 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
469 integer :: m, is, ie, js, je, isd, ied, jsd, jed
470 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
471 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
473 if (.not.
associated(cs))
return
475 if (cs%coupled_tracers)
then
479 call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
480 state%tr_fields, idim=(/isd, is, ie, ied/), &
481 jdim=(/jsd, js, je, jed/) )
493 if (
associated(cs))
then
494 if (
associated(cs%tr))
deallocate(cs%tr)