24 use coupler_types_mod,
only : ind_flux, ind_alpha, ind_csurf
25 use coupler_types_mod,
only : coupler_type_extract_data, coupler_type_set_data
28 implicit none ;
private
30 #include <MOM_memory.h>
32 public register_ocmip2_cfc, initialize_ocmip2_cfc, flux_init_ocmip2_cfc
33 public ocmip2_cfc_column_physics, ocmip2_cfc_surface_state
34 public ocmip2_cfc_stock, ocmip2_cfc_end
37 integer,
parameter :: ntr = 2
41 character(len=200) :: ic_file
44 type(time_type),
pointer :: time => null()
46 real,
pointer,
dimension(:,:,:) :: &
65 real :: cfc11_ic_val = 0.0
66 real :: cfc12_ic_val = 0.0
67 real :: cfc11_land_val = -1.0
68 real :: cfc12_land_val = -1.0
69 logical :: tracers_may_reinit
71 character(len=16) :: cfc11_name
72 character(len=16) :: cfc12_name
74 integer :: ind_cfc_11_flux
76 integer :: ind_cfc_12_flux
92 function register_ocmip2_cfc(HI, GV, param_file, CS, tr_Reg, restart_CS)
105 character(len=40) :: mdl =
"MOM_OCMIP2_CFC"
106 character(len=200) :: inputdir
108 #include "version_variable.h"
109 real,
dimension(:,:,:),
pointer :: tr_ptr => null()
110 real :: a11_dflt(4), a12_dflt(4)
111 real :: d11_dflt(4), d12_dflt(4)
112 real :: e11_dflt(3), e12_dflt(3)
113 character(len=48) :: flux_units
114 logical :: register_ocmip2_cfc
115 integer :: isd, ied, jsd, jed, nz, m
117 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
119 if (
associated(cs))
then
120 call mom_error(warning,
"register_OCMIP2_CFC called with an "// &
121 "associated control structure.")
128 call flux_init_ocmip2_cfc(cs, verbosity=3)
129 if ((cs%ind_cfc_11_flux < 0) .or. (cs%ind_cfc_11_flux < 0))
then
132 call mom_error(warning,
"CFCs are currently only set up to be run in " // &
133 " coupled model configurations, and will be disabled.")
135 register_ocmip2_cfc = .false.
141 call get_param(param_file, mdl,
"CFC_IC_FILE", cs%IC_file, &
142 "The file in which the CFC initial values can be "//&
143 "found, or an empty string for internal initialization.", &
145 if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file,
'/') == 0))
then
147 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
148 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
149 call log_param(param_file, mdl,
"INPUTDIR/CFC_IC_FILE", cs%IC_file)
151 call get_param(param_file, mdl,
"CFC_IC_FILE_IS_Z", cs%Z_IC_file, &
152 "If true, CFC_IC_FILE is in depth space, not layer space", &
154 call get_param(param_file, mdl,
"TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
155 "If true, tracers may go through the initialization code "//&
156 "if they are not found in the restart files. Otherwise "//&
157 "it is a fatal error if tracers are not found in the "//&
158 "restart files of a restarted run.", default=.false.)
162 cs%CFC11_name =
"CFC11" ; cs%CFC12_name =
"CFC12"
163 cs%CFC11_desc = var_desc(cs%CFC11_name,
"mol m-3",
"CFC-11 Concentration", caller=mdl)
164 cs%CFC12_desc = var_desc(cs%CFC12_name,
"mol m-3",
"CFC-12 Concentration", caller=mdl)
165 if (gv%Boussinesq)
then ; flux_units =
"mol s-1"
166 else ; flux_units =
"mol m-3 kg s-1" ;
endif
168 allocate(cs%CFC11(isd:ied,jsd:jed,nz)) ; cs%CFC11(:,:,:) = 0.0
169 allocate(cs%CFC12(isd:ied,jsd:jed,nz)) ; cs%CFC12(:,:,:) = 0.0
175 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
176 tr_desc=cs%CFC11_desc, registry_diags=.true., &
177 flux_units=flux_units, &
178 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
181 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
182 tr_desc=cs%CFC12_desc, registry_diags=.true., &
183 flux_units=flux_units, &
184 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
192 a11_dflt(:) = (/ 3501.8, -210.31, 6.1851, -0.07513 /)
193 a12_dflt(:) = (/ 3845.4, -228.95, 6.1908, -0.06743 /)
194 call get_param(param_file, mdl,
"CFC11_A1", cs%a1_11, &
195 "A coefficient in the Schmidt number of CFC11.", &
196 units=
"nondim", default=a11_dflt(1))
197 call get_param(param_file, mdl,
"CFC11_A2", cs%a2_11, &
198 "A coefficient in the Schmidt number of CFC11.", &
199 units=
"degC-1", default=a11_dflt(2))
200 call get_param(param_file, mdl,
"CFC11_A3", cs%a3_11, &
201 "A coefficient in the Schmidt number of CFC11.", &
202 units=
"degC-2", default=a11_dflt(3))
203 call get_param(param_file, mdl,
"CFC11_A4", cs%a4_11, &
204 "A coefficient in the Schmidt number of CFC11.", &
205 units=
"degC-3", default=a11_dflt(4))
207 call get_param(param_file, mdl,
"CFC12_A1", cs%a1_12, &
208 "A coefficient in the Schmidt number of CFC12.", &
209 units=
"nondim", default=a12_dflt(1))
210 call get_param(param_file, mdl,
"CFC12_A2", cs%a2_12, &
211 "A coefficient in the Schmidt number of CFC12.", &
212 units=
"degC-1", default=a12_dflt(2))
213 call get_param(param_file, mdl,
"CFC12_A3", cs%a3_12, &
214 "A coefficient in the Schmidt number of CFC12.", &
215 units=
"degC-2", default=a12_dflt(3))
216 call get_param(param_file, mdl,
"CFC12_A4", cs%a4_12, &
217 "A coefficient in the Schmidt number of CFC12.", &
218 units=
"degC-3", default=a12_dflt(4))
224 d11_dflt(:) = (/ -229.9261, 319.6552, 119.4471, -1.39165 /)
225 e11_dflt(:) = (/ -0.142382, 0.091459, -0.0157274 /)
226 d12_dflt(:) = (/ -218.0971, 298.9702, 113.8049, -1.39165 /)
227 e12_dflt(:) = (/ -0.143566, 0.091015, -0.0153924 /)
229 call get_param(param_file, mdl,
"CFC11_D1", cs%d1_11, &
230 "A coefficient in the solubility of CFC11.", &
231 units=
"none", default=d11_dflt(1))
232 call get_param(param_file, mdl,
"CFC11_D2", cs%d2_11, &
233 "A coefficient in the solubility of CFC11.", &
234 units=
"hK", default=d11_dflt(2))
235 call get_param(param_file, mdl,
"CFC11_D3", cs%d3_11, &
236 "A coefficient in the solubility of CFC11.", &
237 units=
"none", default=d11_dflt(3))
238 call get_param(param_file, mdl,
"CFC11_D4", cs%d4_11, &
239 "A coefficient in the solubility of CFC11.", &
240 units=
"hK-2", default=d11_dflt(4))
241 call get_param(param_file, mdl,
"CFC11_E1", cs%e1_11, &
242 "A coefficient in the solubility of CFC11.", &
243 units=
"PSU-1", default=e11_dflt(1))
244 call get_param(param_file, mdl,
"CFC11_E2", cs%e2_11, &
245 "A coefficient in the solubility of CFC11.", &
246 units=
"PSU-1 hK-1", default=e11_dflt(2))
247 call get_param(param_file, mdl,
"CFC11_E3", cs%e3_11, &
248 "A coefficient in the solubility of CFC11.", &
249 units=
"PSU-1 hK-2", default=e11_dflt(3))
251 call get_param(param_file, mdl,
"CFC12_D1", cs%d1_12, &
252 "A coefficient in the solubility of CFC12.", &
253 units=
"none", default=d12_dflt(1))
254 call get_param(param_file, mdl,
"CFC12_D2", cs%d2_12, &
255 "A coefficient in the solubility of CFC12.", &
256 units=
"hK", default=d12_dflt(2))
257 call get_param(param_file, mdl,
"CFC12_D3", cs%d3_12, &
258 "A coefficient in the solubility of CFC12.", &
259 units=
"none", default=d12_dflt(3))
260 call get_param(param_file, mdl,
"CFC12_D4", cs%d4_12, &
261 "A coefficient in the solubility of CFC12.", &
262 units=
"hK-2", default=d12_dflt(4))
263 call get_param(param_file, mdl,
"CFC12_E1", cs%e1_12, &
264 "A coefficient in the solubility of CFC12.", &
265 units=
"PSU-1", default=e12_dflt(1))
266 call get_param(param_file, mdl,
"CFC12_E2", cs%e2_12, &
267 "A coefficient in the solubility of CFC12.", &
268 units=
"PSU-1 hK-1", default=e12_dflt(2))
269 call get_param(param_file, mdl,
"CFC12_E3", cs%e3_12, &
270 "A coefficient in the solubility of CFC12.", &
271 units=
"PSU-1 hK-2", default=e12_dflt(3))
274 cs%restart_CSp => restart_cs
276 register_ocmip2_cfc = .true.
277 end function register_ocmip2_cfc
281 subroutine flux_init_ocmip2_cfc(CS, verbosity)
285 integer,
optional,
intent(in) :: verbosity
288 character(len=128) :: default_ice_restart_file =
'ice_ocmip2_cfc.res.nc'
289 character(len=128) :: default_ocean_restart_file =
'ocmip2_cfc.res.nc'
290 integer :: ind_flux(2)
294 ind_flux(1) = aof_set_coupler_flux(
'cfc_11_flux', &
295 flux_type =
'air_sea_gas_flux', implementation =
'ocmip2', &
296 param = (/ 9.36e-07, 9.7561e-06 /), &
297 ice_restart_file = default_ice_restart_file, &
298 ocean_restart_file = default_ocean_restart_file, &
299 caller =
"register_OCMIP2_CFC", verbosity=verbosity)
300 ind_flux(2) = aof_set_coupler_flux(
'cfc_12_flux', &
301 flux_type =
'air_sea_gas_flux', implementation =
'ocmip2', &
302 param = (/ 9.36e-07, 9.7561e-06 /), &
303 ice_restart_file = default_ice_restart_file, &
304 ocean_restart_file = default_ocean_restart_file, &
305 caller =
"register_OCMIP2_CFC", verbosity=verbosity)
307 if (
present(cs))
then ;
if (
associated(cs))
then
308 cs%ind_cfc_11_flux = ind_flux(1)
309 cs%ind_cfc_12_flux = ind_flux(2)
312 end subroutine flux_init_ocmip2_cfc
315 subroutine initialize_ocmip2_cfc(restart, day, G, GV, US, h, diag, OBC, CS, &
317 logical,
intent(in) :: restart
319 type(time_type),
target,
intent(in) :: day
323 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
325 type(
diag_ctrl),
target,
intent(in) :: diag
338 logical :: from_file = .false.
340 if (.not.
associated(cs))
return
345 if (.not.restart .or. (cs%tracers_may_reinit .and. &
347 call init_tracer_cfc(h, cs%CFC11, cs%CFC11_name, cs%CFC11_land_val, &
348 cs%CFC11_IC_val, g, us, cs)
350 if (.not.restart .or. (cs%tracers_may_reinit .and. &
352 call init_tracer_cfc(h, cs%CFC12, cs%CFC12_name, cs%CFC12_land_val, &
353 cs%CFC12_IC_val, g, us, cs)
355 if (
associated(obc))
then
359 end subroutine initialize_ocmip2_cfc
362 subroutine init_tracer_cfc(h, tr, name, land_val, IC_val, G, US, CS)
365 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
366 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: tr
367 character(len=*),
intent(in) :: name
368 real,
intent(in) :: land_val
369 real,
intent(in) :: IC_val
376 integer :: i, j, k, is, ie, js, je, nz
377 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
379 if (len_trim(cs%IC_file) > 0)
then
382 call mom_error(fatal,
"initialize_OCMIP2_CFC: Unable to open "//cs%IC_file)
383 if (cs%Z_IC_file)
then
384 ok = tracer_z_init(tr, h, cs%IC_file, name, g, us)
386 ok = tracer_z_init(tr, h, cs%IC_file, trim(name), g, us)
387 if (.not.ok)
call mom_error(fatal,
"initialize_OCMIP2_CFC: "//&
388 "Unable to read "//trim(name)//
" from "//&
389 trim(cs%IC_file)//
".")
395 do k=1,nz ;
do j=js,je ;
do i=is,ie
396 if (g%mask2dT(i,j) < 0.5)
then
401 enddo ;
enddo ;
enddo
404 end subroutine init_tracer_cfc
409 subroutine ocmip2_cfc_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
410 evap_CFL_limit, minimum_forcing_depth)
413 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
415 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
417 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
421 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
425 type(
forcing),
intent(in) :: fluxes
427 real,
intent(in) :: dt
430 real,
optional,
intent(in) :: evap_cfl_limit
432 real,
optional,
intent(in) :: minimum_forcing_depth
444 real :: c1(szi_(g),szk_(g))
445 real,
dimension(SZI_(G),SZJ_(G)) :: &
448 real,
pointer,
dimension(:,:,:) :: cfc11 => null(), cfc12 => null()
449 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
450 integer :: i, j, k, m, is, ie, js, je, nz, idim(4), jdim(4)
452 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
453 idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)
455 if (.not.
associated(cs))
return
457 cfc11 => cs%CFC11 ; cfc12 => cs%CFC12
462 call coupler_type_extract_data(fluxes%tr_fluxes, cs%ind_cfc_11_flux, ind_flux, &
463 cfc11_flux, -gv%Rho0, idim=idim, jdim=jdim)
464 call coupler_type_extract_data(fluxes%tr_fluxes, cs%ind_cfc_12_flux, ind_flux, &
465 cfc12_flux, -gv%Rho0, idim=idim, jdim=jdim)
469 if (
present(evap_cfl_limit) .and.
present(minimum_forcing_depth))
then
470 do k=1,nz ;
do j=js,je ;
do i=is,ie
471 h_work(i,j,k) = h_old(i,j,k)
472 enddo ;
enddo ;
enddo
473 call applytracerboundaryfluxesinout(g, gv, cfc11, dt, fluxes, h_work, &
474 evap_cfl_limit, minimum_forcing_depth)
475 call tracer_vertdiff(h_work, ea, eb, dt, cfc11, g, gv, sfc_flux=cfc11_flux)
477 do k=1,nz ;
do j=js,je ;
do i=is,ie
478 h_work(i,j,k) = h_old(i,j,k)
479 enddo ;
enddo ;
enddo
480 call applytracerboundaryfluxesinout(g, gv, cfc12, dt, fluxes, h_work, &
481 evap_cfl_limit, minimum_forcing_depth)
482 call tracer_vertdiff(h_work, ea, eb, dt, cfc12, g, gv, sfc_flux=cfc12_flux)
484 call tracer_vertdiff(h_old, ea, eb, dt, cfc11, g, gv, sfc_flux=cfc11_flux)
485 call tracer_vertdiff(h_old, ea, eb, dt, cfc12, g, gv, sfc_flux=cfc12_flux)
490 end subroutine ocmip2_cfc_column_physics
495 function ocmip2_cfc_stock(h, stocks, G, GV, CS, names, units, stock_index)
498 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
500 real,
dimension(:),
intent(out) :: stocks
504 character(len=*),
dimension(:),
intent(out) :: names
505 character(len=*),
dimension(:),
intent(out) :: units
506 integer,
optional,
intent(in) :: stock_index
508 integer :: ocmip2_cfc_stock
512 integer :: i, j, k, is, ie, js, je, nz
513 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
516 if (.not.
associated(cs))
return
518 if (
present(stock_index))
then ;
if (stock_index > 0)
then
525 call query_vardesc(cs%CFC11_desc, name=names(1), units=units(1), caller=
"OCMIP2_CFC_stock")
526 call query_vardesc(cs%CFC12_desc, name=names(2), units=units(2), caller=
"OCMIP2_CFC_stock")
527 units(1) = trim(units(1))//
" kg" ; units(2) = trim(units(2))//
" kg"
529 stocks(1) = 0.0 ; stocks(2) = 0.0
530 do k=1,nz ;
do j=js,je ;
do i=is,ie
531 mass = g%mask2dT(i,j) * g%areaT(i,j) * h(i,j,k)
532 stocks(1) = stocks(1) + cs%CFC11(i,j,k) * mass
533 stocks(2) = stocks(2) + cs%CFC12(i,j,k) * mass
534 enddo ;
enddo ;
enddo
535 stocks(1) = gv%H_to_kg_m2 * stocks(1)
536 stocks(2) = gv%H_to_kg_m2 * stocks(2)
540 end function ocmip2_cfc_stock
544 subroutine ocmip2_cfc_surface_state(state, h, G, CS)
546 type(
surface),
intent(inout) :: state
548 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
554 real,
dimension(SZI_(G),SZJ_(G)) :: &
566 integer :: i, j, m, is, ie, js, je, idim(4), jdim(4)
568 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
569 idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)
571 if (.not.
associated(cs))
return
573 do j=js,je ;
do i=is,ie
574 ta = max(0.01, (state%SST(i,j) + 273.15) * 0.01)
575 sal = state%SSS(i,j) ; sst = state%SST(i,j)
580 alpha_11 = exp(cs%d1_11 + cs%d2_11/ta + cs%d3_11*log(ta) + cs%d4_11*ta**2 +&
581 sal * ((cs%e3_11 * ta + cs%e2_11) * ta + cs%e1_11)) * &
582 1.0e-09 * g%mask2dT(i,j)
583 alpha_12 = exp(cs%d1_12 + cs%d2_12/ta + cs%d3_12*log(ta) + cs%d4_12*ta**2 +&
584 sal * ((cs%e3_12 * ta + cs%e2_12) * ta + cs%e1_12)) * &
585 1.0e-09 * g%mask2dT(i,j)
588 sc_11 = cs%a1_11 + sst * (cs%a2_11 + sst * (cs%a3_11 + sst * cs%a4_11)) * &
590 sc_12 = cs%a1_12 + sst * (cs%a2_12 + sst * (cs%a3_12 + sst * cs%a4_12)) * &
593 sc_no_term = sqrt(660.0 / (abs(sc_11) + 1.0e-30))
594 cfc11_alpha(i,j) = alpha_11 * sc_no_term
595 cfc11_csurf(i,j) = cs%CFC11(i,j,1) * sc_no_term
597 sc_no_term = sqrt(660.0 / (abs(sc_12) + 1.0e-30))
598 cfc12_alpha(i,j) = alpha_12 * sc_no_term
599 cfc12_csurf(i,j) = cs%CFC12(i,j,1) * sc_no_term
604 call coupler_type_set_data(cfc11_alpha, cs%ind_cfc_11_flux, ind_alpha, &
605 state%tr_fields, idim=idim, jdim=jdim)
606 call coupler_type_set_data(cfc11_csurf, cs%ind_cfc_11_flux, ind_csurf, &
607 state%tr_fields, idim=idim, jdim=jdim)
608 call coupler_type_set_data(cfc12_alpha, cs%ind_cfc_12_flux, ind_alpha, &
609 state%tr_fields, idim=idim, jdim=jdim)
610 call coupler_type_set_data(cfc12_csurf, cs%ind_cfc_12_flux, ind_csurf, &
611 state%tr_fields, idim=idim, jdim=jdim)
613 end subroutine ocmip2_cfc_surface_state
616 subroutine ocmip2_cfc_end(CS)
624 if (
associated(cs))
then
625 if (
associated(cs%CFC11))
deallocate(cs%CFC11)
626 if (
associated(cs%CFC12))
deallocate(cs%CFC12)
630 end subroutine ocmip2_cfc_end