MOM_OCMIP2_CFC.F90
1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Simulates CFCs using the OCMIP2 protocols
6module mom_ocmip2_cfc
7
9use mom_coupler_types, only : extract_coupler_type_data, set_coupler_type_data
10use mom_coupler_types, only : atmos_ocn_coupler_flux
11use mom_diag_mediator, only : diag_ctrl
12use mom_error_handler, only : mom_error, fatal, warning
13use mom_file_parser, only : get_param, log_param, log_version, param_file_type
14use mom_forcing_type, only : forcing
15use mom_hor_index, only : hor_index_type
16use mom_grid, only : ocean_grid_type
17use mom_io, only : file_exists, mom_read_data, slasher
18use mom_io, only : vardesc, var_desc, query_vardesc
19use mom_open_boundary, only : ocean_obc_type
20use mom_restart, only : query_initialized, set_initialized, mom_restart_cs
21use mom_spatial_means, only : global_mass_int_efp
22use mom_sponge, only : set_up_sponge_field, sponge_cs
23use mom_time_manager, only : time_type
24use mom_tracer_registry, only : register_tracer, tracer_registry_type
25use mom_tracer_diabatic, only : tracer_vertdiff, applytracerboundaryfluxesinout
26use mom_tracer_z_init, only : tracer_z_init
27use mom_unit_scaling, only : unit_scale_type
28use mom_variables, only : surface
29use mom_verticalgrid, only : verticalgrid_type
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
37public ocmip2_cfc_stock, ocmip2_cfc_end
38
39!> The control structure for the OCMPI2_CFC tracer package
40type, public :: ocmip2_cfc_cs ; private
41 character(len=200) :: ic_file !< The file in which the CFC initial values can
42 !! be found, or an empty string for internal initilaization.
43 logical :: z_ic_file !< If true, the IC_file is in Z-space. The default is false..
44 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
45 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM6 tracer registry
46 real, pointer, dimension(:,:,:) :: &
47 cfc11 => null(), & !< The CFC11 concentration [mol m-3].
48 cfc12 => null() !< The CFC12 concentration [mol m-3].
49 ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12.
50 !>@{ Coefficients used in the CFC11 and CFC12 solubility calculation
51 real :: a1_11, a1_12 ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [nondim]
52 real :: a2_11, a2_12 ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-1]
53 real :: a3_11, a3_12 ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-2]
54 real :: a4_11, a4_12 ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-3]
55
56 real :: d1_11, d1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [nondim]
57 real :: d2_11, d2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-1]
58 real :: d3_11, d3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [log(hectoKelvin)-1]
59 real :: d4_11, d4_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-2]
60
61 real :: e1_11, e1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1]
62 real :: e2_11, e2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1 hectoKelvin-1]
63 real :: e3_11, e3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-2 hectoKelvin-2]
64 !>@}
65 real :: cfc11_ic_val = 0.0 !< The initial value assigned to CFC11 [mol m-3].
66 real :: cfc12_ic_val = 0.0 !< The initial value assigned to CFC12 [mol m-3].
67 real :: cfc11_land_val = -1.0 !< The value of CFC11 used where land is masked out [mol m-3].
68 real :: cfc12_land_val = -1.0 !< The value of CFC12 used where land is masked out [mol m-3].
69 logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code
70 !! if they are not found in the restart files.
71 character(len=16) :: cfc11_name !< CFC11 variable name
72 character(len=16) :: cfc12_name !< CFC12 variable name
73
74 integer :: ind_cfc_11_flux !< Index returned by atmos_ocn_coupler_flux that is used to
75 !! pack and unpack surface boundary condition arrays.
76 integer :: ind_cfc_12_flux !< Index returned by atmos_ocn_coupler_flux that is used to
77 !! pack and unpack surface boundary condition arrays.
78
79 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate
80 !! the timing of diagnostic output.
81 type(mom_restart_cs), pointer :: restart_csp => null() !< Model restart control structure
82
83 ! The following vardesc types contain a package of metadata about each tracer.
84 type(vardesc) :: cfc11_desc !< A set of metadata for the CFC11 tracer
85 type(vardesc) :: cfc12_desc !< A set of metadata for the CFC12 tracer
86end type ocmip2_cfc_cs
87
88contains
89
90!> Register the OCMIP2 CFC tracers to be used with MOM and read the parameters
91!! that are used with this tracer package
92function register_ocmip2_cfc(HI, GV, param_file, CS, tr_Reg, restart_CS)
93 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
94 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
95 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
96 type(ocmip2_cfc_cs), pointer :: cs !< A pointer that is set to point to the control
97 !! structure for this module.
98 type(tracer_registry_type), &
99 pointer :: tr_reg !< A pointer to the tracer registry.
100 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
101
102 ! Local variables
103 character(len=40) :: mdl = "MOM_OCMIP2_CFC" ! This module's name.
104 character(len=200) :: inputdir ! The directory where NetCDF input files are.
105 ! This include declares and sets the variable "version".
106# include "version_variable.h"
107 real, dimension(:,:,:), pointer :: tr_ptr => null() ! A pointer to a CFC tracer [mol m-3]
108 real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients
109 real :: d11_dflt(4), d12_dflt(4) ! in the expressions for the solubility and
110 real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers [various units by element].
111 character(len=48) :: flux_units ! The units for tracer fluxes.
112 logical :: register_ocmip2_cfc
113 integer :: isd, ied, jsd, jed, nz
114
115 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
116
117 if (associated(cs)) then
118 call mom_error(fatal, "register_OCMIP2_CFC called with an "// &
119 "associated control structure.")
120 endif
121 allocate(cs)
122
123 ! This call sets default properties for the air-sea CFC fluxes and obtains the
124 ! indicies for the CFC11 and CFC12 flux coupling.
125 call flux_init_ocmip2_cfc(cs, verbosity=3)
126 if ((cs%ind_cfc_11_flux < 0) .or. (cs%ind_cfc_12_flux < 0)) then
127 ! This is most likely to happen with the dummy version of atmos_ocn_coupler_flux
128 ! used in ocean-only runs.
129 call mom_error(warning, "CFCs are currently only set up to be run in " // &
130 " coupled model configurations, and will be disabled.")
131 deallocate(cs)
132 register_ocmip2_cfc = .false.
133 return
134 endif
135
136 ! Read all relevant parameters and write them to the model log.
137 call log_version(param_file, mdl, version, "")
138 call get_param(param_file, mdl, "CFC_IC_FILE", cs%IC_file, &
139 "The file in which the CFC initial values can be "//&
140 "found, or an empty string for internal initialization.", &
141 default=" ")
142 if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file,'/') == 0)) then
143 ! Add the directory if CS%IC_file is not already a complete path.
144 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
145 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
146 call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", cs%IC_file)
147 endif
148 call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", cs%Z_IC_file, &
149 "If true, CFC_IC_FILE is in depth space, not layer space", &
150 default=.false.)
151 call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
152 "If true, tracers may go through the initialization code "//&
153 "if they are not found in the restart files. Otherwise "//&
154 "it is a fatal error if tracers are not found in the "//&
155 "restart files of a restarted run.", default=.false.)
156
157 ! The following vardesc types contain a package of metadata about each tracer,
158 ! including, the name; units; longname; and grid information.
159 cs%CFC11_name = "CFC11" ; cs%CFC12_name = "CFC12"
160 cs%CFC11_desc = var_desc(cs%CFC11_name,"mol m-3","CFC-11 Concentration", caller=mdl)
161 cs%CFC12_desc = var_desc(cs%CFC12_name,"mol m-3","CFC-12 Concentration", caller=mdl)
162 if (gv%Boussinesq) then ; flux_units = "mol s-1"
163 else ; flux_units = "mol m-3 kg s-1" ; endif
164
165 allocate(cs%CFC11(isd:ied,jsd:jed,nz), source=0.0)
166 allocate(cs%CFC12(isd:ied,jsd:jed,nz), source=0.0)
167
168 ! This pointer assignment is needed to force the compiler not to do a copy in
169 ! the registration calls. Curses on the designers and implementers of F90.
170 tr_ptr => cs%CFC11
171 ! Register CFC11 for horizontal advection, diffusion, and restarts.
172 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
173 tr_desc=cs%CFC11_desc, registry_diags=.true., &
174 flux_units=flux_units, &
175 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
176 ! Do the same for CFC12
177 tr_ptr => cs%CFC12
178 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
179 tr_desc=cs%CFC12_desc, registry_diags=.true., &
180 flux_units=flux_units, &
181 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
182
183 ! Set and read the various empirical coefficients.
184
185!-----------------------------------------------------------------------
186! Default Schmidt number coefficients for CFC11 (_11) and CFC12 (_12) are given
187! by Zheng et al (1998), JGR vol 103, C1.
188!-----------------------------------------------------------------------
189 a11_dflt(:) = (/ 3501.8, -210.31, 6.1851, -0.07513 /)
190 a12_dflt(:) = (/ 3845.4, -228.95, 6.1908, -0.06743 /)
191 call get_param(param_file, mdl, "CFC11_A1", cs%a1_11, &
192 "A coefficient in the Schmidt number of CFC11.", &
193 units="nondim", default=a11_dflt(1))
194 call get_param(param_file, mdl, "CFC11_A2", cs%a2_11, &
195 "A coefficient in the Schmidt number of CFC11.", &
196 units="degC-1", default=a11_dflt(2))
197 call get_param(param_file, mdl, "CFC11_A3", cs%a3_11, &
198 "A coefficient in the Schmidt number of CFC11.", &
199 units="degC-2", default=a11_dflt(3))
200 call get_param(param_file, mdl, "CFC11_A4", cs%a4_11, &
201 "A coefficient in the Schmidt number of CFC11.", &
202 units="degC-3", default=a11_dflt(4))
203
204 call get_param(param_file, mdl, "CFC12_A1", cs%a1_12, &
205 "A coefficient in the Schmidt number of CFC12.", &
206 units="nondim", default=a12_dflt(1))
207 call get_param(param_file, mdl, "CFC12_A2", cs%a2_12, &
208 "A coefficient in the Schmidt number of CFC12.", &
209 units="degC-1", default=a12_dflt(2))
210 call get_param(param_file, mdl, "CFC12_A3", cs%a3_12, &
211 "A coefficient in the Schmidt number of CFC12.", &
212 units="degC-2", default=a12_dflt(3))
213 call get_param(param_file, mdl, "CFC12_A4", cs%a4_12, &
214 "A coefficient in the Schmidt number of CFC12.", &
215 units="degC-3", default=a12_dflt(4))
216
217!-----------------------------------------------------------------------
218! Solubility coefficients for alpha in mol/l/atm for CFC11 (_11) and CFC12 (_12)
219! after Warner and Weiss (1985) DSR, vol 32.
220!-----------------------------------------------------------------------
221 d11_dflt(:) = (/ -229.9261, 319.6552, 119.4471, -1.39165 /)
222 e11_dflt(:) = (/ -0.142382, 0.091459, -0.0157274 /)
223 d12_dflt(:) = (/ -218.0971, 298.9702, 113.8049, -1.39165 /)
224 e12_dflt(:) = (/ -0.143566, 0.091015, -0.0153924 /)
225
226 call get_param(param_file, mdl, "CFC11_D1", cs%d1_11, &
227 "A coefficient in the solubility of CFC11.", &
228 units="none", default=d11_dflt(1))
229 call get_param(param_file, mdl, "CFC11_D2", cs%d2_11, &
230 "A coefficient in the solubility of CFC11.", &
231 units="hK", default=d11_dflt(2))
232 call get_param(param_file, mdl, "CFC11_D3", cs%d3_11, &
233 "A coefficient in the solubility of CFC11.", &
234 units="none", default=d11_dflt(3))
235 call get_param(param_file, mdl, "CFC11_D4", cs%d4_11, &
236 "A coefficient in the solubility of CFC11.", &
237 units="hK-2", default=d11_dflt(4))
238 call get_param(param_file, mdl, "CFC11_E1", cs%e1_11, &
239 "A coefficient in the solubility of CFC11.", &
240 units="PSU-1", default=e11_dflt(1))
241 call get_param(param_file, mdl, "CFC11_E2", cs%e2_11, &
242 "A coefficient in the solubility of CFC11.", &
243 units="PSU-1 hK-1", default=e11_dflt(2))
244 call get_param(param_file, mdl, "CFC11_E3", cs%e3_11, &
245 "A coefficient in the solubility of CFC11.", &
246 units="PSU-1 hK-2", default=e11_dflt(3))
247
248 call get_param(param_file, mdl, "CFC12_D1", cs%d1_12, &
249 "A coefficient in the solubility of CFC12.", &
250 units="none", default=d12_dflt(1))
251 call get_param(param_file, mdl, "CFC12_D2", cs%d2_12, &
252 "A coefficient in the solubility of CFC12.", &
253 units="hK", default=d12_dflt(2))
254 call get_param(param_file, mdl, "CFC12_D3", cs%d3_12, &
255 "A coefficient in the solubility of CFC12.", &
256 units="none", default=d12_dflt(3))
257 call get_param(param_file, mdl, "CFC12_D4", cs%d4_12, &
258 "A coefficient in the solubility of CFC12.", &
259 units="hK-2", default=d12_dflt(4))
260 call get_param(param_file, mdl, "CFC12_E1", cs%e1_12, &
261 "A coefficient in the solubility of CFC12.", &
262 units="PSU-1", default=e12_dflt(1))
263 call get_param(param_file, mdl, "CFC12_E2", cs%e2_12, &
264 "A coefficient in the solubility of CFC12.", &
265 units="PSU-1 hK-1", default=e12_dflt(2))
266 call get_param(param_file, mdl, "CFC12_E3", cs%e3_12, &
267 "A coefficient in the solubility of CFC12.", &
268 units="PSU-1 hK-2", default=e12_dflt(3))
269
270 cs%tr_Reg => tr_reg
271 cs%restart_CSp => restart_cs
272
273 register_ocmip2_cfc = .true.
274end function register_ocmip2_cfc
275
276!> This subroutine initializes the air-sea CFC fluxes, and optionally returns
277!! the indicies of these fluxes. It can safely be called multiple times.
278subroutine flux_init_ocmip2_cfc(CS, verbosity)
279 type(ocmip2_cfc_cs), optional, pointer :: cs !< An optional pointer to the control structure
280 !! for this module; if not present, the flux indicies
281 !! are not stored.
282 integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity.
283
284 ! These can be overridden later in via the field manager?
285 character(len=128) :: default_ice_restart_file = 'ice_ocmip2_cfc.res.nc'
286 character(len=128) :: default_ocean_restart_file = 'ocmip2_cfc.res.nc'
287 integer :: ind_flux(2) ! Integer indices of the fluxes
288
289 ! These calls obtain the indices for the CFC11 and CFC12 flux coupling. They
290 ! can safely be called multiple times.
291 ind_flux(1) = atmos_ocn_coupler_flux('cfc_11_flux', &
292 flux_type='air_sea_gas_flux', implementation='ocmip2', &
293 param=(/ 9.36e-07, 9.7561e-06 /), &
294 ice_restart_file=default_ice_restart_file, &
295 ocean_restart_file=default_ocean_restart_file, &
296 caller="register_OCMIP2_CFC", verbosity=verbosity)
297 ind_flux(2) = atmos_ocn_coupler_flux('cfc_12_flux', &
298 flux_type='air_sea_gas_flux', implementation='ocmip2', &
299 param=(/ 9.36e-07, 9.7561e-06 /), &
300 ice_restart_file=default_ice_restart_file, &
301 ocean_restart_file=default_ocean_restart_file, &
302 caller="register_OCMIP2_CFC", verbosity=verbosity)
303
304 if (present(cs)) then ; if (associated(cs)) then
305 cs%ind_cfc_11_flux = ind_flux(1)
306 cs%ind_cfc_12_flux = ind_flux(2)
307 endif ; endif
308
309end subroutine flux_init_ocmip2_cfc
310
311!> Initialize the OCMP2 CFC tracer fields and set up the tracer output.
312subroutine initialize_ocmip2_cfc(restart, day, G, GV, US, h, diag, OBC, CS, &
313 sponge_CSp)
314 logical, intent(in) :: restart !< .true. if the fields have already been
315 !! read from a restart file.
316 type(time_type), target, intent(in) :: day !< Time of the start of the run.
317 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
318 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
319 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
320 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
321 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
322 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
323 !! diagnostic output.
324 type(ocean_obc_type), pointer :: obc !< This open boundary condition type
325 !! specifies whether, where, and what
326 !! open boundary conditions are used.
327 type(ocmip2_cfc_cs), pointer :: cs !< The control structure returned by a
328 !! previous call to register_OCMIP2_CFC.
329 type(sponge_cs), pointer :: sponge_csp !< A pointer to the control structure for
330 !! the sponges, if they are in use.
331 !! Otherwise this may be unassociated.
332
333 if (.not.associated(cs)) return
334
335 cs%Time => day
336 cs%diag => diag
337
338 if (.not.restart .or. (cs%tracers_may_reinit .and. &
339 .not.query_initialized(cs%CFC11, cs%CFC11_name, cs%restart_CSp))) then
340 call init_tracer_cfc(h, cs%CFC11, cs%CFC11_name, cs%CFC11_land_val, &
341 cs%CFC11_IC_val, g, gv, us, cs)
342 call set_initialized(cs%CFC11, cs%CFC11_name, cs%restart_CSp)
343 endif
344
345 if (.not.restart .or. (cs%tracers_may_reinit .and. &
346 .not.query_initialized(cs%CFC12, cs%CFC12_name, cs%restart_CSp))) then
347 call init_tracer_cfc(h, cs%CFC12, cs%CFC12_name, cs%CFC12_land_val, &
348 cs%CFC12_IC_val, g, gv, us, cs)
349 call set_initialized(cs%CFC12, cs%CFC12_name, cs%restart_CSp)
350 endif
351
352 if (associated(obc)) then
353 ! Steal from updated DOME in the fullness of time.
354 endif
355
356end subroutine initialize_ocmip2_cfc
357
358!>This subroutine initializes a tracer array.
359subroutine init_tracer_cfc(h, tr, name, land_val, IC_val, G, GV, US, CS)
360 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
361 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
362 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
363 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
364 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: tr !< The CFC tracer concentration array [mol m-3]
365 character(len=*), intent(in) :: name !< The tracer name
366 real, intent(in) :: land_val !< A value the tracer takes over land [mol m-3]
367 real, intent(in) :: IC_val !< The initial condition value for
368 !! the CRC tracer [mol m-3]
369 type(ocmip2_cfc_cs), pointer :: CS !< The control structure returned by a
370 !! previous call to register_OCMIP2_CFC.
371
372 ! This subroutine initializes a tracer array.
373
374 logical :: OK
375 integer :: i, j, k, is, ie, js, je, nz
376 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
377
378 if (len_trim(cs%IC_file) > 0) then
379 ! Read the tracer concentrations from a netcdf file.
380 if (.not.file_exists(cs%IC_file, g%Domain)) &
381 call mom_error(fatal, "initialize_OCMIP2_CFC: Unable to open "//cs%IC_file)
382 if (cs%Z_IC_file) then
383 ok = tracer_z_init(tr, h, cs%IC_file, name, g, gv, us)
384 if (.not.ok) then
385 ok = tracer_z_init(tr, h, cs%IC_file, trim(name), g, gv, us)
386 if (.not.ok) call mom_error(fatal,"initialize_OCMIP2_CFC: "//&
387 "Unable to read "//trim(name)//" from "//&
388 trim(cs%IC_file)//".")
389 endif
390 else
391 call mom_read_data(cs%IC_file, trim(name), tr, g%Domain)
392 endif
393 else
394 do k=1,nz ; do j=js,je ; do i=is,ie
395 if (g%mask2dT(i,j) < 0.5) then
396 tr(i,j,k) = land_val
397 else
398 tr(i,j,k) = ic_val
399 endif
400 enddo ; enddo ; enddo
401 endif
402
403end subroutine init_tracer_cfc
404
405!> This subroutine applies diapycnal diffusion, souces and sinks and any other column
406!! tracer physics or chemistry to the OCMIP2 CFC tracers.
407!! CFCs are relatively simple, as they are passive tracers with only a surface flux as a source.
408subroutine ocmip2_cfc_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
409 evap_CFL_limit, minimum_forcing_depth)
410 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
411 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
412 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
413 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
414 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
415 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
416 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
417 intent(in) :: ea !< an array to which the amount of fluid entrained
418 !! from the layer above during this call will be
419 !! added [H ~> m or kg m-2].
420 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
421 intent(in) :: eb !< an array to which the amount of fluid entrained
422 !! from the layer below during this call will be
423 !! added [H ~> m or kg m-2].
424 type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic
425 !! and tracer forcing fields. Unused fields have NULL ptrs.
426 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
427 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
428 type(ocmip2_cfc_cs), pointer :: cs !< The control structure returned by a
429 !! previous call to register_OCMIP2_CFC.
430 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
431 !! be fluxed out of the top layer in a timestep [nondim]
432 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
433 !! fluxes can be applied [H ~> m or kg m-2]
434! This subroutine applies diapycnal diffusion and any other column
435! tracer physics or chemistry to the tracers from this file.
436! CFCs are relatively simple, as they are passive tracers. with only a surface
437! flux as a source.
438
439! The arguments to this subroutine are redundant in that
440! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
441
442 ! Local variables
443 real, dimension(SZI_(G),SZJ_(G)) :: &
444 cfc11_flux, & ! The fluxes of CFC11 and CFC12 into the ocean, in unscaled units of
445 cfc12_flux ! CFC concentrations times a vertical mass flux [mol R Z m-3 T-1 ~> mol kg m-3 s-1]
446 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
447 integer :: i, j, k, is, ie, js, je, nz, idim(4), jdim(4)
448
449 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
450 idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)
451
452 if (.not.associated(cs)) return
453
454 ! These two calls unpack the fluxes from the input arrays.
455 ! The -GV%Rho0 changes the sign convention of the flux and with the scaling factors changes
456 ! the units of the flux from [conc m s-1] to [conc R Z T-1 ~> conc kg m-2 s-1].
457 call extract_coupler_type_data(fluxes%tr_fluxes, cs%ind_cfc_11_flux, cfc11_flux, &
458 scale_factor=-gv%Rho0*us%m_to_Z*us%T_to_s, idim=idim, jdim=jdim, turns=g%HI%turns)
459 call extract_coupler_type_data(fluxes%tr_fluxes, cs%ind_cfc_12_flux, cfc12_flux, &
460 scale_factor=-gv%Rho0*us%m_to_Z*us%T_to_s, idim=idim, jdim=jdim, turns=g%HI%turns)
461
462 ! Use a tridiagonal solver to determine the concentrations after the
463 ! surface source is applied and diapycnal advection and diffusion occurs.
464 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
465 do k=1,nz ;do j=js,je ; do i=is,ie
466 h_work(i,j,k) = h_old(i,j,k)
467 enddo ; enddo ; enddo
468 call applytracerboundaryfluxesinout(g, gv, cs%CFC11, dt, fluxes, h_work, &
469 evap_cfl_limit, minimum_forcing_depth)
470 call tracer_vertdiff(h_work, ea, eb, dt, cs%CFC11, g, gv, sfc_flux=cfc11_flux)
471
472 do k=1,nz ;do j=js,je ; do i=is,ie
473 h_work(i,j,k) = h_old(i,j,k)
474 enddo ; enddo ; enddo
475 call applytracerboundaryfluxesinout(g, gv, cs%CFC12, dt, fluxes, h_work, &
476 evap_cfl_limit, minimum_forcing_depth)
477 call tracer_vertdiff(h_work, ea, eb, dt, cs%CFC12, g, gv, sfc_flux=cfc12_flux)
478 else
479 call tracer_vertdiff(h_old, ea, eb, dt, cs%CFC11, g, gv, sfc_flux=cfc11_flux)
480 call tracer_vertdiff(h_old, ea, eb, dt, cs%CFC12, g, gv, sfc_flux=cfc12_flux)
481 endif
482
483 ! Write out any desired diagnostics from tracer sources & sinks here.
484
485end subroutine ocmip2_cfc_column_physics
486
487!> This function calculates the mass-weighted integral of all tracer stocks,
488!! returning the number of stocks it has calculated. If the stock_index
489!! is present, only the stock corresponding to that coded index is returned.
490function ocmip2_cfc_stock(h, stocks, G, GV, CS, names, units, stock_index)
491 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
492 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
493 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
494 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
495 type(efp_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each
496 !! tracer, in kg times concentration units [kg conc]
497 type(ocmip2_cfc_cs), pointer :: cs !< The control structure returned by a
498 !! previous call to register_OCMIP2_CFC.
499 character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated.
500 character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated.
501 integer, optional, intent(in) :: stock_index !< The coded index of a specific
502 !! stock being sought.
503 integer :: ocmip2_cfc_stock !< The number of stocks calculated here.
504
505
506 ocmip2_cfc_stock = 0
507 if (.not.associated(cs)) return
508
509 if (present(stock_index)) then ; if (stock_index > 0) then
510 ! Check whether this stock is available from this routine.
511
512 ! No stocks from this routine are being checked yet. Return 0.
513 return
514 endif ; endif
515
516 call query_vardesc(cs%CFC11_desc, name=names(1), units=units(1), caller="OCMIP2_CFC_stock")
517 call query_vardesc(cs%CFC12_desc, name=names(2), units=units(2), caller="OCMIP2_CFC_stock")
518 units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"
519
520 stocks(1) = global_mass_int_efp(h, g, gv, cs%CFC11, on_pe_only=.true.)
521 stocks(2) = global_mass_int_efp(h, g, gv, cs%CFC12, on_pe_only=.true.)
522
523 ocmip2_cfc_stock = 2
524
525end function ocmip2_cfc_stock
526
527!> This subroutine extracts the surface CFC concentrations and other fields that
528!! are shared with the atmosphere to calculate CFC fluxes.
529subroutine ocmip2_cfc_surface_state(sfc_state, h, G, GV, US, CS)
530 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
531 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
532 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
533 !! describe the surface state of the ocean.
534 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
535 intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
536 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
537 type(ocmip2_cfc_cs), pointer :: cs !< The control structure returned by a previous
538 !! call to register_OCMIP2_CFC.
539
540 ! Local variables
541 real, dimension(SZI_(G),SZJ_(G)) :: &
542 cfc11_csurf, & ! The CFC-11 surface concentrations times the Schmidt number term [mol m-3].
543 cfc12_csurf, & ! The CFC-12 surface concentrations times the Schmidt number term [mol m-3].
544 cfc11_alpha, & ! The CFC-11 solubility [mol m-3 pptv-1].
545 cfc12_alpha ! The CFC-12 solubility [mol m-3 pptv-1].
546 real :: ta ! Absolute sea surface temperature [hectoKelvin] (Why use such bizzare units?)
547 real :: sal ! Surface salinity [PSU].
548 real :: sst ! Sea surface temperature [degC].
549 real :: alpha_11 ! The solubility of CFC 11 [mol m-3 pptv-1].
550 real :: alpha_12 ! The solubility of CFC 12 [mol m-3 pptv-1].
551 real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12 [nondim].
552 real :: sc_no_term ! A term related to the Schmidt number [nondim].
553 integer :: i, j, is, ie, js, je, idim(4), jdim(4)
554
555 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
556 idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)
557
558 if (.not.associated(cs)) return
559
560 do j=js,je ; do i=is,ie
561 ta = max(0.01, (us%C_to_degC*sfc_state%SST(i,j) + 273.15) * 0.01) ! Why is this in hectoKelvin?
562 sal = us%S_to_ppt*sfc_state%SSS(i,j) ; sst = us%C_to_degC*sfc_state%SST(i,j)
563 ! Calculate solubilities using Warner and Weiss (1985) DSR, vol 32.
564 ! The final result is in mol/cm3/pptv (1 part per trillion 1e-12)
565 ! Use Bullister and Wisegavger for CCl4.
566 ! The factor 1.e-09 converts from mol/(l * atm) to mol/(m3 * pptv).
567 alpha_11 = exp(cs%d1_11 + cs%d2_11/ta + cs%d3_11*log(ta) + cs%d4_11*ta**2 +&
568 sal * ((cs%e3_11 * ta + cs%e2_11) * ta + cs%e1_11)) * &
569 1.0e-09 * g%mask2dT(i,j)
570 alpha_12 = exp(cs%d1_12 + cs%d2_12/ta + cs%d3_12*log(ta) + cs%d4_12*ta**2 +&
571 sal * ((cs%e3_12 * ta + cs%e2_12) * ta + cs%e1_12)) * &
572 1.0e-09 * g%mask2dT(i,j)
573 ! Calculate Schmidt numbers using coefficients given by
574 ! Zheng et al (1998), JGR vol 103, C1.
575 sc_11 = cs%a1_11 + sst * (cs%a2_11 + sst * (cs%a3_11 + sst * cs%a4_11)) * &
576 g%mask2dT(i,j)
577 sc_12 = cs%a1_12 + sst * (cs%a2_12 + sst * (cs%a3_12 + sst * cs%a4_12)) * &
578 g%mask2dT(i,j)
579 ! The abs here is to avoid NaNs. The model should be failing at this point.
580 sc_no_term = sqrt(660.0 / (abs(sc_11) + 1.0e-30))
581 cfc11_alpha(i,j) = alpha_11 * sc_no_term
582 cfc11_csurf(i,j) = cs%CFC11(i,j,1) * sc_no_term
583
584 sc_no_term = sqrt(660.0 / (abs(sc_12) + 1.0e-30))
585 cfc12_alpha(i,j) = alpha_12 * sc_no_term
586 cfc12_csurf(i,j) = cs%CFC12(i,j,1) * sc_no_term
587 enddo ; enddo
588
589 ! These calls load these values into the appropriate arrays in the
590 ! coupler-type structure.
591 call set_coupler_type_data(cfc11_alpha, cs%ind_cfc_11_flux, sfc_state%tr_fields, &
592 solubility=.true., idim=idim, jdim=jdim, turns=g%HI%turns)
593 call set_coupler_type_data(cfc11_csurf, cs%ind_cfc_11_flux, sfc_state%tr_fields, &
594 idim=idim, jdim=jdim, turns=g%HI%turns)
595 call set_coupler_type_data(cfc12_alpha, cs%ind_cfc_12_flux, sfc_state%tr_fields, &
596 solubility=.true., idim=idim, jdim=jdim, turns=g%HI%turns)
597 call set_coupler_type_data(cfc12_csurf, cs%ind_cfc_12_flux, sfc_state%tr_fields, &
598 idim=idim, jdim=jdim, turns=g%HI%turns)
599
600end subroutine ocmip2_cfc_surface_state
601
602!> Deallocate any memory associated with the OCMIP2 CFC tracer package
603subroutine ocmip2_cfc_end(CS)
604 type(ocmip2_cfc_cs), pointer :: cs !< The control structure returned by a
605 !! previous call to register_OCMIP2_CFC.
606! This subroutine deallocates the memory owned by this module.
607! Argument: CS - The control structure returned by a previous call to
608! register_OCMIP2_CFC.
609
610 if (associated(cs)) then
611 if (associated(cs%CFC11)) deallocate(cs%CFC11)
612 if (associated(cs%CFC12)) deallocate(cs%CFC12)
613
614 deallocate(cs)
615 endif
616end subroutine ocmip2_cfc_end
617
618
619!> \namespace mom_ocmip2_cfc
620!!
621!! By Robert Hallberg, 2007
622!!
623!! This module contains the code that is needed to set
624!! up and use CFC-11 and CFC-12 in a fully coupled or ice-ocean model
625!! context using the OCMIP2 protocols
626
627end module mom_ocmip2_cfc