MOM_CFC_cap.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 atmospheric pressure, wind speed and sea ice cover
6!! provided via cap (only NUOPC cap is implemented so far).
7module mom_cfc_cap
8
9use mom_coms, only : efp_type
10use mom_debugging, only : hchksum
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
16use mom_grid, only : ocean_grid_type
18use mom_io, only : file_exists, mom_read_data, slasher
19use mom_io, only : vardesc, var_desc, query_vardesc, stdout
20use mom_tracer_registry, only : tracer_type
22use mom_restart, only : query_initialized, set_initialized, mom_restart_cs
24use mom_time_manager, only : time_type, increment_date
25use mom_interpolate, only : external_field, init_external_field, time_interp_external
33
34implicit none ; private
35
36#include <MOM_memory.h>
37
41
42integer, parameter :: ntr = 2 !< the number of tracers in this module.
43
44!> Contains the concentration array, surface flux, a pointer to Tr in Tr_reg,
45!! and some metadata for a single CFC tracer
46type, private :: cfc_tracer_data
47 type(vardesc) :: desc !< A set of metadata for the tracer
48 real :: ic_val = 0.0 !< The initial value assigned to the tracer [mol kg-1].
49 real :: land_val = -1.0 !< The value of the tracer used where land is
50 !! masked out [mol kg-1].
51 character(len=32) :: name !< Tracer variable name
52 integer :: id_cmor = -1 !< Diagnostic id
53 integer :: id_sfc_flux = -1 !< Surface flux id
54 real, pointer, dimension(:,:,:) :: conc !< The tracer concentration [mol kg-1].
55 real, pointer, dimension(:,:) :: sfc_flux !< Surface flux [CU R Z T-1 ~> mol m-2 s-1]
56 type(tracer_type), pointer :: tr_ptr !< pointer to tracer inside Tr_reg
57end type cfc_tracer_data
58
59!> The control structure for the CFC_cap tracer package
60type, public :: cfc_cap_cs ; private
61 logical :: debug !< If true, write verbose checksums for debugging purposes.
62 character(len=200) :: ic_file !< The file in which the CFC initial values can
63 !! be found, or an empty string for internal initilaization.
64 logical :: z_ic_file !< If true, the IC_file is in Z-space. The default is false.
65 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
66 type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM6 tracer registry
67 logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code
68 !! if they are not found in the restart files.
69 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate
70 !! the timing of diagnostic output.
71 type(mom_restart_cs), pointer :: restart_csp => null() !< Model restart control structure
72
73 type(cfc_tracer_data), dimension(NTR) :: cfc_data !< per-tracer parameters / metadata
74 integer :: cfc_bc_year_offset = 0 !< offset to add to model time to get time value used in CFC_BC_file
75 type(external_field) :: cfc11_atm_nh_handle !< Handle for time-interpolated CFC11 atm NH
76 type(external_field) :: cfc11_atm_sh_handle !< Handle for time-interpolated CFC11 atm SH
77 type(external_field) :: cfc12_atm_nh_handle !< Handle for time-interpolated CFC12 atm NH
78 type(external_field) :: cfc12_atm_sh_handle !< Handle for time-interpolated CFC12 atm SH
79end type cfc_cap_cs
80
81contains
82
83!> Register the CFCs to be used with MOM and read the parameters
84!! that are used with this tracer package
85function register_cfc_cap(HI, GV, param_file, CS, tr_Reg, restart_CS)
86 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
87 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
88 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
89 type(cfc_cap_cs), pointer :: cs !< A pointer that is set to point to the control
90 !! structure for this module.
92 pointer :: tr_reg !< A pointer to the tracer registry.
93 type(mom_restart_cs), target, intent(inout) :: restart_cs !< MOM restart control struct
94
95 ! Local variables
96 character(len=40) :: mdl = "MOM_CFC_cap" ! This module's name.
97 ! This include declares and sets the variable "version".
98# include "version_variable.h"
99 character(len=200) :: inputdir ! The directory where NetCDF input files are.
100 real, dimension(:,:,:), pointer :: tr_ptr => null() ! A pointer to a CFC tracer [mol kg-1]
101 character(len=200) :: cfc_bc_file ! filename with cfc11 and cfc12 data
102 character(len=30) :: cfc_bc_var_name ! varname of field in CFC_BC_file
103 character :: m2char
104 logical :: register_cfc_cap
105 integer :: isd, ied, jsd, jed, nz, m
106 integer :: cfc_bc_data_year ! specific year in CFC BC data calendar
107 integer :: cfc_bc_model_year ! model year corresponding to CFC_BC_data_year
108
109 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
110
111 if (associated(cs)) then
112 call mom_error(fatal, "register_CFC_cap called with an "// &
113 "associated control structure.")
114 endif
115 allocate(cs)
116
117 ! Read all relevant parameters and write them to the model log.
118 call log_version(param_file, mdl, version, "")
119 call get_param(param_file, mdl, "DEBUG", cs%debug, &
120 "If true, write out verbose debugging data.", &
121 default=.false., debuggingparam=.true.)
122 call get_param(param_file, mdl, "CFC_IC_FILE", cs%IC_file, &
123 "The file in which the CFC initial values can be "//&
124 "found, or an empty string for internal initialization.", &
125 default=" ")
126 if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file, '/') == 0)) then
127 ! Add the directory if CS%IC_file is not already a complete path.
128 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
129 cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
130 call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", cs%IC_file, &
131 "full path of CFC_IC_FILE")
132 endif
133 call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", cs%Z_IC_file, &
134 "If true, CFC_IC_FILE is in depth space, not layer space", &
135 default=.false.)
136 call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
137 "If true, tracers may go through the initialization code "//&
138 "if they are not found in the restart files. Otherwise "//&
139 "it is a fatal error if tracers are not found in the "//&
140 "restart files of a restarted run.", default=.false.)
141 do m=1,ntr
142 write(m2char, "(I1)") m
143 call get_param(param_file, mdl, "CFC1"//m2char//"_IC_VAL", cs%CFC_data(m)%IC_val, &
144 "Value that CFC_1"//m2char//" is set to when it is not read from a file.", &
145 units="mol kg-1", default=0.0)
146 enddo
147
148 ! the following params are not used in this module. Instead, they are used in
149 ! the cap but are logged here to keep all the CFC cap params together.
150 call get_param(param_file, mdl, "CFC_BC_FILE", cfc_bc_file, &
151 "The file in which the CFC-11 and CFC-12 atm concentrations can be "//&
152 "found (units must be parts per trillion).", default=" ")
153 if (len_trim(cfc_bc_file) == 0) then
154 call mom_error(fatal, "CFC_BC_FILE must be specified if USE_CFC_CAP=.true.")
155 endif
156 if (scan(cfc_bc_file, '/') == 0) then
157 ! Add the directory if CFC_BC_file is not already a complete path.
158 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
159 cfc_bc_file = trim(slasher(inputdir))//trim(cfc_bc_file)
160 call log_param(param_file, mdl, "INPUTDIR/CFC_BC_FILE", cfc_bc_file, &
161 "full path of CFC_BC_FILE")
162 endif
163
164 call get_param(param_file, mdl, "CFC_BC_DATA_YEAR", cfc_bc_data_year, &
165 "Specific year in CFC_BC_FILE data calendar", default=2000)
166 call get_param(param_file, mdl, "CFC_BC_MODEL_YEAR", cfc_bc_model_year, &
167 "Model year corresponding to CFC_BC_MODEL_YEAR", default=2000)
168 cs%CFC_BC_year_offset = cfc_bc_data_year - cfc_bc_model_year
169
170 call get_param(param_file, mdl, "CFC11_NH_VARIABLE", cfc_bc_var_name, &
171 "Variable name of NH CFC-11 atm mole fraction in CFC_BC_FILE.", &
172 default="cfc11_nh")
173 cs%cfc11_atm_nh_handle = init_external_field(cfc_bc_file, cfc_bc_var_name)
174
175 call get_param(param_file, mdl, "CFC11_SH_VARIABLE", cfc_bc_var_name, &
176 "Variable name of SH CFC-11 atm mole fraction in CFC_BC_FILE.", &
177 default="cfc11_sh")
178 cs%cfc11_atm_sh_handle = init_external_field(cfc_bc_file, cfc_bc_var_name)
179
180 call get_param(param_file, mdl, "CFC12_NH_VARIABLE", cfc_bc_var_name, &
181 "Variable name of NH CFC-12 atm mole fraction in CFC_BC_FILE.", &
182 default="cfc12_nh")
183 cs%cfc12_atm_nh_handle = init_external_field(cfc_bc_file, cfc_bc_var_name)
184
185 call get_param(param_file, mdl, "CFC12_SH_VARIABLE", cfc_bc_var_name, &
186 "Variable name of SH CFC-12 atm mole fraction in CFC_BC_FILE.", &
187 default="cfc12_sh")
188 cs%cfc12_atm_sh_handle = init_external_field(cfc_bc_file, cfc_bc_var_name)
189! domain=G%Domain%mpp_domain)
190
191 ! The following vardesc types contain a package of metadata about each tracer,
192 ! including, the name; units; longname; and grid information.
193 do m=1,ntr
194 write(m2char, "(I1)") m
195 write(cs%CFC_data(m)%name, "(2A)") "CFC_1", m2char
196 cs%CFC_data(m)%desc = var_desc(cs%CFC_data(m)%name, &
197 "mol kg-1", &
198 "Moles Per Unit Mass of CFC-1"//m2char//" in sea water", &
199 caller=mdl)
200
201 allocate(cs%CFC_data(m)%conc(isd:ied,jsd:jed,nz), source=0.0)
202 allocate(cs%CFC_data(m)%sfc_flux(isd:ied,jsd:jed), source=0.0)
203
204 ! This pointer assignment is needed to force the compiler not to do a copy in
205 ! the registration calls. Curses on the designers and implementers of F90.
206 tr_ptr => cs%CFC_data(m)%conc
207 ! Register CFC tracer for horizontal advection, diffusion, and restarts.
208 call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
209 tr_desc=cs%CFC_data(m)%desc, registry_diags=.true., &
210 restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit, &
211 tr_out=cs%CFC_data(m)%tr_ptr)
212 enddo
213
214 cs%tr_Reg => tr_reg
215 cs%restart_CSp => restart_cs
216 register_cfc_cap = .true.
217
218end function register_cfc_cap
219
220!> Initialize the CFC tracer fields and set up the tracer output.
221subroutine initialize_cfc_cap(restart, day, G, GV, US, h, diag, OBC, CS)
222 logical, intent(in) :: restart !< .true. if the fields have already been
223 !! read from a restart file.
224 type(time_type), target, intent(in) :: day !< Time of the start of the run.
225 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
226 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
227 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
228 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
229 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
230 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
231 !! diagnostic output.
232 type(ocean_obc_type), pointer :: obc !< This open boundary condition type
233 !! specifies whether, where, and what
234 !! open boundary conditions are used.
235 type(cfc_cap_cs), pointer :: cs !< The control structure returned by a
236 !! previous call to register_CFC_cap.
237
238 ! local variables
239 integer :: m
240 character :: m2char
241
242 if (.not.associated(cs)) return
243
244 cs%Time => day
245 cs%diag => diag
246
247 do m=1,ntr
248 if (.not.restart .or. (cs%tracers_may_reinit .and. &
249 .not.query_initialized(cs%CFC_data(m)%conc, cs%CFC_data(m)%name, cs%restart_CSp))) then
250 call init_tracer_cfc(h, cs%CFC_data(m)%conc, cs%CFC_data(m)%name, cs%CFC_data(m)%land_val, &
251 cs%CFC_data(m)%IC_val, g, gv, us, cs)
252 call set_initialized(cs%CFC_data(m)%conc, cs%CFC_data(m)%name, cs%restart_CSp)
253 endif
254
255 ! cmor diagnostics
256 ! units for cfc11_flux and cfc12_flux are [Conc R Z T-1 ~> mol m-2 s-1]
257 ! CFC11 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/42625c97b8fe75124a345962c4430982.html
258 ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html
259 ! CFC12 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/3ab8e10027d7014f18f9391890369235.html
260 ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html
261 write(m2char, "(I1)") m
262 cs%CFC_data(m)%id_cmor = register_diag_field('ocean_model', &
263 'cfc1'//m2char, diag%axesTL, day, &
264 'Mole Concentration of CFC1'//m2char//' in Sea Water', 'mol m-3', &
265 conversion=gv%Rho0*us%R_to_kg_m3)
266
267 cs%CFC_data(m)%id_sfc_flux = register_diag_field('ocean_model', &
268 'cfc1'//m2char//'_flux', diag%axesT1, day, &
269 'Gas exchange flux of CFC1'//m2char//' into the ocean ', &
270 'mol m-2 s-1', conversion=us%RZ_T_to_kg_m2s, &
271 cmor_field_name='fgcfc1'//m2char, &
272 cmor_long_name='Surface Downward CFC1'//m2char//' Flux', &
273 cmor_standard_name='surface_downward_cfc1'//m2char//'_flux')
274 enddo
275
276
277 if (associated(obc)) then
278 ! Steal from updated DOME in the fullness of time.
279 ! GMM: TODO this must be coded if we intend to use this module in regional applications
280 endif
281
282end subroutine initialize_cfc_cap
283
284!>This subroutine initializes a tracer array.
285subroutine init_tracer_cfc(h, tr, name, land_val, IC_val, G, GV, US, CS)
286 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
287 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
288 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
289 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
290 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: tr !< The tracer concentration array [mol kg-1]
291 character(len=*), intent(in) :: name !< The tracer name
292 real, intent(in) :: land_val !< A value the tracer takes over land [mol kg-1]
293 real, intent(in) :: IC_val !< The initial condition value for the
294 !! tracer [mol kg-1]
295 type(cfc_cap_cs), pointer :: CS !< The control structure returned by a
296 !! previous call to register_CFC_cap.
297
298 ! local variables
299 logical :: OK
300 integer :: i, j, k, is, ie, js, je, nz
301 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
302
303 if (len_trim(cs%IC_file) > 0) then
304 ! Read the tracer concentrations from a netcdf file.
305 if (.not.file_exists(cs%IC_file, g%Domain)) &
306 call mom_error(fatal, "initialize_CFC_cap: Unable to open "//cs%IC_file)
307 if (cs%Z_IC_file) then
308 ok = tracer_z_init(tr, h, cs%IC_file, name, g, gv, us)
309 if (.not.ok) then
310 ok = tracer_z_init(tr, h, cs%IC_file, trim(name), g, gv, us)
311 if (.not.ok) call mom_error(fatal,"initialize_CFC_cap: "//&
312 "Unable to read "//trim(name)//" from "//&
313 trim(cs%IC_file)//".")
314 endif
315 else
316 call mom_read_data(cs%IC_file, trim(name), tr, g%Domain)
317 endif
318 else
319 do k=1,nz ; do j=js,je ; do i=is,ie
320 if (g%mask2dT(i,j) < 0.5) then
321 tr(i,j,k) = land_val
322 else
323 tr(i,j,k) = ic_val
324 endif
325 enddo ; enddo ; enddo
326 endif
327
328end subroutine init_tracer_cfc
329
330!> Applies diapycnal diffusion, souces and sinks and any other column
331!! tracer physics to the CFC cap tracers. CFCs are relatively simple,
332!! as they are passive tracers with only a surface flux as a source.
333subroutine cfc_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, KPP_CSp, &
334 nonLocalTrans, evap_CFL_limit, minimum_forcing_depth)
335 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
336 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
337 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
338 intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
339 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
340 intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
341 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
342 intent(in) :: ea !< an array to which the amount of fluid entrained
343 !! from the layer above during this call will be
344 !! added [H ~> m or kg m-2].
345 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
346 intent(in) :: eb !< an array to which the amount of fluid entrained
347 !! from the layer below during this call will be
348 !! added [H ~> m or kg m-2].
349 type(forcing), intent(in) :: fluxes!< A structure containing pointers to thermodynamic
350 !! and tracer forcing fields. Unused fields have NULL ptrs.
351 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]
352 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
353 type(cfc_cap_cs), pointer :: cs !< The control structure returned by a
354 !! previous call to register_CFC_cap.
355 type(kpp_cs), optional, pointer :: kpp_csp !< KPP control structure
356 real, optional, intent(in) :: nonlocaltrans(:,:,:) !< Non-local transport [nondim]
357 real, optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
358 !! be fluxed out of the top layer in a timestep [nondim]
359 real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
360 !! fluxes can be applied [H ~> m or kg m-2]
361
362 ! The arguments to this subroutine are redundant in that
363 ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)
364
365 ! Local variables
366 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
367 integer :: i, j, k, is, ie, js, je, nz, m
368
369 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
370
371 if (.not.associated(cs)) return
372
373 ! Compute KPP nonlocal term if necessary
374 if (present(kpp_csp)) then
375 if (associated(kpp_csp) .and. present(nonlocaltrans)) then
376 do m=1,ntr
377 call kpp_nonlocaltransport(kpp_csp, g, gv, h_old, nonlocaltrans, &
378 cs%CFC_data(m)%sfc_flux(:,:), dt, cs%diag, &
379 cs%CFC_data(m)%tr_ptr, cs%CFC_data(m)%conc(:,:,:), &
380 flux_scale=gv%RZ_to_H)
381 enddo
382 endif
383 endif
384
385 ! Use a tridiagonal solver to determine the concentrations after the
386 ! surface source is applied and diapycnal advection and diffusion occurs.
387 if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
388 do m=1,ntr
389 do k=1,nz ;do j=js,je ; do i=is,ie
390 h_work(i,j,k) = h_old(i,j,k)
391 enddo ; enddo ; enddo
392 call applytracerboundaryfluxesinout(g, gv, cs%CFC_data(m)%conc, dt, fluxes, h_work, &
393 evap_cfl_limit, minimum_forcing_depth)
394 call tracer_vertdiff(h_work, ea, eb, dt, cs%CFC_data(m)%conc, g, gv, &
395 sfc_flux=cs%CFC_data(m)%sfc_flux)
396 enddo
397 else
398 do m=1,ntr
399 call tracer_vertdiff(h_old, ea, eb, dt, cs%CFC_data(m)%conc, g, gv, &
400 sfc_flux=cs%CFC_data(m)%sfc_flux)
401 enddo
402 endif
403
404 ! If needed, write out any desired diagnostics from tracer sources & sinks here.
405 do m=1,ntr
406 if (cs%CFC_data(m)%id_cmor > 0) &
407 call post_data(cs%CFC_data(m)%id_cmor, cs%CFC_data(m)%conc, cs%diag)
408
409 if (cs%CFC_data(m)%id_sfc_flux > 0) &
410 call post_data(cs%CFC_data(m)%id_sfc_flux, cs%CFC_data(m)%sfc_flux, cs%diag)
411 enddo
412
413end subroutine cfc_cap_column_physics
414
415
416!> Calculates the mass-weighted integral of all tracer stocks,
417!! returning the number of stocks it has calculated. If the stock_index
418!! is present, only the stock corresponding to that coded index is returned.
419function cfc_cap_stock(h, stocks, G, GV, CS, names, units, stock_index)
420 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
421 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
422 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
423 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
424 type(efp_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each
425 !! tracer, in kg times concentration units [kg conc]
426 type(cfc_cap_cs), pointer :: cs !< The control structure returned by a
427 !! previous call to register_CFC_cap.
428 character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated.
429 character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated.
430 integer, optional, intent(in) :: stock_index !< The coded index of a specific
431 !! stock being sought.
432 integer :: cfc_cap_stock !< The number of stocks calculated here.
433
434 ! Local variables
435 real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1]
436 real :: mass ! The cell volume or mass [H L2 ~> m3 or kg]
437 integer :: i, j, k, is, ie, js, je, nz, m
438 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
439
440 cfc_cap_stock = 0
441 if (.not.associated(cs)) return
442
443 if (present(stock_index)) then ; if (stock_index > 0) then
444 ! Check whether this stock is available from this routine.
445
446 ! No stocks from this routine are being checked yet. Return 0.
447 return
448 endif ; endif
449
450 do m=1,ntr
451 call query_vardesc(cs%CFC_data(m)%desc, name=names(m), units=units(m), caller="CFC_cap_stock")
452 units(m) = trim(units(m))//" kg"
453 stocks(m) = global_mass_int_efp(h, g, gv, cs%CFC_data(m)%conc, on_pe_only=.true.)
454 enddo
455
456 cfc_cap_stock = ntr
457
458end function cfc_cap_stock
459
460!> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM
461!! concentration, and calculating the solubility, Schmidt number, and gas exchange.
462subroutine cfc_cap_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, Rho0, CS)
463 type(surface), intent(in ) :: sfc_state !< A structure containing fields
464 !! that describe the surface state of the ocean.
465 type(forcing), intent(inout) :: fluxes !< A structure containing pointers
466 !! to thermodynamic and tracer forcing fields. Unused fields
467 !! have NULL ptrs.
468 type(time_type), intent(in) :: day_start !< Start time of the fluxes.
469 type(time_type), intent(in) :: day_interval !< Length of time over which these
470 !! fluxes will be applied.
471 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
472 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
473 real, intent(in) :: rho0 !< The mean ocean density [R ~> kg m-3]
474 type(cfc_cap_cs), pointer :: cs !< The control structure returned by a
475 !! previous call to register_CFC_cap.
476
477 ! Local variables
478 type(time_type) :: time_external ! time value used in CFC_BC_file
479 real, dimension(SZI_(G),SZJ_(G)) :: &
480 kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1].
481 kw, & ! gas transfer velocity [Z T-1 ~> m s-1].
482 cair, & ! The surface gas concentration in equilibrium with the atmosphere
483 ! (saturation concentration) [mol kg-1].
484 cfc11_atm, & ! CFC11 atm mole fraction [pico mol/mol]
485 cfc12_atm ! CFC12 atm mole fraction [pico mol/mol]
486 real :: cfc11_atm_nh ! NH value for cfc11_atm [pico mol/mol]
487 real :: cfc11_atm_sh ! SH value for cfc11_atm [pico mol/mol]
488 real :: cfc12_atm_nh ! NH value for cfc12_atm [pico mol/mol]
489 real :: cfc12_atm_sh ! SH value for cfc12_atm [pico mol/mol]
490 real :: ta ! Absolute sea surface temperature [hectoKelvin]
491 real :: sal ! Surface salinity [PSU].
492 real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1].
493 real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1].
494 real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12 [nondim].
495 real :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2] = [Z T L-2 ~> s m-1]
496 real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm [atm Pa-1].
497 real :: press_to_atm ! converts from model pressure units to atm [atm T2 R-1 L-2 ~> atm Pa-1]
498 integer :: i, j, is, ie, js, je, m
499
500 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
501
502 ! Time_external = increment_date(day_start + day_interval/2, years=CS%CFC_BC_year_offset)
503 time_external = increment_date(day_start, years=cs%CFC_BC_year_offset)
504
505 ! CFC11 atm mole fraction, convert from ppt (pico mol/mol) to mol/mol
506 call time_interp_external(cs%cfc11_atm_nh_handle, time_external, cfc11_atm_nh)
507 cfc11_atm_nh = cfc11_atm_nh * 1.0e-12
508 call time_interp_external(cs%cfc11_atm_sh_handle, time_external, cfc11_atm_sh)
509 cfc11_atm_sh = cfc11_atm_sh * 1.0e-12
510
511 ! CFC12 atm mole fraction, convert from ppt (pico mol/mol) to mol/mol
512 call time_interp_external(cs%cfc12_atm_nh_handle, time_external, cfc12_atm_nh)
513 cfc12_atm_nh = cfc12_atm_nh * 1.0e-12
514 call time_interp_external(cs%cfc12_atm_sh_handle, time_external, cfc12_atm_sh)
515 cfc12_atm_sh = cfc12_atm_sh * 1.0e-12
516
517 !---------------------------------------------------------------------
518 ! Gas exchange/piston velocity parameter
519 !---------------------------------------------------------------------
520 ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014
521 ! = 6.97e-7 [m/s s^2/m^2] [Z T-1 T2 L-2] = [Z T L-2 ~> s m-1]
522 kw_coeff = (us%m_to_Z*us%s_to_T*us%L_to_m**2) * 6.97e-7
523
524 ! set unit conversion factors
525 press_to_atm = us%R_to_kg_m3*us%L_T_to_m_s**2 * pa_to_atm
526
527 do j=js,je ; do i=is,ie
528 if (g%geoLatT(i,j) < -10.0) then
529 cfc11_atm(i,j) = cfc11_atm_sh
530 cfc12_atm(i,j) = cfc12_atm_sh
531 elseif (g%geoLatT(i,j) <= 10.0) then
532 cfc11_atm(i,j) = cfc11_atm_sh + &
533 (0.05 * g%geoLatT(i,j) + 0.5) * (cfc11_atm_nh - cfc11_atm_sh)
534 cfc12_atm(i,j) = cfc12_atm_sh + &
535 (0.05 * g%geoLatT(i,j) + 0.5) * (cfc12_atm_nh - cfc12_atm_sh)
536 else
537 cfc11_atm(i,j) = cfc11_atm_nh
538 cfc12_atm(i,j) = cfc12_atm_nh
539 endif
540 enddo ; enddo
541
542 do j=js,je ; do i=is,ie
543 ! ta in hectoKelvin
544 ta = max(0.01, (us%C_to_degC*sfc_state%SST(i,j) + 273.15) * 0.01)
545 sal = us%S_to_ppt*sfc_state%SSS(i,j)
546
547 ! Calculate solubilities
548 call get_solubility(alpha_11, alpha_12, ta, sal , g%mask2dT(i,j))
549
550 ! Calculate Schmidt numbers using coefficients given by
551 ! Wanninkhof (2014); doi:10.4319/lom.2014.12.351.
552 call comp_cfc_schmidt(us%C_to_degC*sfc_state%SST(i,j), sc_11, sc_12)
553
554 kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j))
555
556 ! air concentrations and cfcs BC's fluxes
557 ! CFC flux units: [mol kg-1 R Z T-1 ~> mol m-2 s-1]
558 kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_11)
559 cair(i,j) = press_to_atm * alpha_11 * cfc11_atm(i,j) * fluxes%p_surf_full(i,j)
560 cs%CFC_data(1)%sfc_flux(i,j) = kw(i,j) * (cair(i,j) - cs%CFC_data(1)%conc(i,j,1)) * rho0
561
562 kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_12)
563 cair(i,j) = press_to_atm * alpha_12 * cfc12_atm(i,j) * fluxes%p_surf_full(i,j)
564 cs%CFC_data(2)%sfc_flux(i,j) = kw(i,j) * (cair(i,j) - cs%CFC_data(2)%conc(i,j,1)) * rho0
565 enddo ; enddo
566
567 if (cs%debug) then
568 do m=1,ntr
569 call hchksum(cs%CFC_data(m)%sfc_flux, trim(cs%CFC_data(m)%name)//" sfc_flux", g%HI, &
570 unscale=us%RZ_T_to_kg_m2s)
571 enddo
572 endif
573
574end subroutine cfc_cap_set_forcing
575
576!> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32.
577subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask)
578 real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1]
579 real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1]
580 real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin]
581 real, intent(in ) :: sal !< Surface salinity [PSU].
582 real, intent(in ) :: mask !< ocean mask [nondim]
583
584 ! Local variables
585
586 ! Coefficients for calculating CFC11 solubilities
587 ! from Table 5 in Warner and Weiss (1985) DSR, vol 32.
588
589 real, parameter :: d1_11 = -232.0411 ! [nondim]
590 real, parameter :: d2_11 = 322.5546 ! [hectoKelvin-1]
591 real, parameter :: d3_11 = 120.4956 ! [log(hectoKelvin)-1]
592 real, parameter :: d4_11 = -1.39165 ! [hectoKelvin-2]
593
594 real, parameter :: e1_11 = -0.146531 ! [PSU-1]
595 real, parameter :: e2_11 = 0.093621 ! [PSU-1 hectoKelvin-1]
596 real, parameter :: e3_11 = -0.0160693 ! [PSU-2 hectoKelvin-2]
597
598 ! Coefficients for calculating CFC12 solubilities
599 ! from Table 5 in Warner and Weiss (1985) DSR, vol 32.
600
601 real, parameter :: d1_12 = -220.2120 ! [nondim]
602 real, parameter :: d2_12 = 301.8695 ! [hectoKelvin-1]
603 real, parameter :: d3_12 = 114.8533 ! [log(hectoKelvin)-1]
604 real, parameter :: d4_12 = -1.39165 ! [hectoKelvin-2]
605
606 real, parameter :: e1_12 = -0.147718 ! [PSU-1]
607 real, parameter :: e2_12 = 0.093175 ! [PSU-1 hectoKelvin-1]
608 real, parameter :: e3_12 = -0.0157340 ! [PSU-2 hectoKelvin-2]
609
610 real :: factor ! introduce units to result [mol kg-1 atm-1]
611
612 ! Eq. 9 from Warner and Weiss (1985) DSR, vol 32.
613 factor = 1.0
614 alpha_11 = exp(d1_11 + d2_11/ta + d3_11*log(ta) + d4_11*ta**2 +&
615 sal * ((e3_11 * ta + e2_11) * ta + e1_11)) * &
616 factor * mask
617 alpha_12 = exp(d1_12 + d2_12/ta + d3_12*log(ta) + d4_12*ta**2 +&
618 sal * ((e3_12 * ta + e2_12) * ta + e1_12)) * &
619 factor * mask
620
621end subroutine get_solubility
622
623
624!> Compute Schmidt numbers of CFCs following Wanninkhof (2014); doi:10.4319/lom.2014.12.351
625!! Range of validity of fit is -2:40.
626subroutine comp_cfc_schmidt(sst_in, cfc11_sc, cfc12_sc)
627 real, intent(in) :: sst_in !< The sea surface temperature [degC].
628 real, intent(inout) :: cfc11_sc !< Schmidt number of CFC11 [nondim].
629 real, intent(inout) :: cfc12_sc !< Schmidt number of CFC12 [nondim].
630
631 !local variables
632 real , parameter :: a_11 = 3579.2 ! CFC11 Schmidt number fit coefficient [nondim]
633 real , parameter :: b_11 = -222.63 ! CFC11 Schmidt number fit coefficient [degC-1]
634 real , parameter :: c_11 = 7.5749 ! CFC11 Schmidt number fit coefficient [degC-2]
635 real , parameter :: d_11 = -0.14595 ! CFC11 Schmidt number fit coefficient [degC-3]
636 real , parameter :: e_11 = 0.0011874 ! CFC11 Schmidt number fit coefficient [degC-4]
637 real , parameter :: a_12 = 3828.1 ! CFC12 Schmidt number fit coefficient [nondim]
638 real , parameter :: b_12 = -249.86 ! CFC12 Schmidt number fit coefficient [degC-1]
639 real , parameter :: c_12 = 8.7603 ! CFC12 Schmidt number fit coefficient [degC-2]
640 real , parameter :: d_12 = -0.1716 ! CFC12 Schmidt number fit coefficient [degC-3]
641 real , parameter :: e_12 = 0.001408 ! CFC12 Schmidt number fit coefficient [degC-4]
642 real :: sst ! A range-limited sea surface temperature [degC]
643
644
645 ! clip SST to avoid bad values
646 sst = max(-2.0, min(40.0, sst_in))
647 cfc11_sc = a_11 + sst * (b_11 + sst * (c_11 + sst * (d_11 + sst * e_11)))
648 cfc12_sc = a_12 + sst * (b_12 + sst * (c_12 + sst * (d_12 + sst * e_12)))
649
650end subroutine comp_cfc_schmidt
651
652!> Deallocate any memory associated with the CFC cap tracer package
653subroutine cfc_cap_end(CS)
654 type(cfc_cap_cs), pointer :: cs !< The control structure returned by a
655 !! previous call to register_CFC_cap.
656
657 ! local variables
658 integer :: m
659
660 if (associated(cs)) then
661 do m=1,ntr
662 if (associated(cs%CFC_data(m)%conc)) deallocate(cs%CFC_data(m)%conc)
663 if (associated(cs%CFC_data(m)%sfc_flux)) deallocate(cs%CFC_data(m)%sfc_flux)
664 enddo
665
666 deallocate(cs)
667 endif
668end subroutine cfc_cap_end
669
670!> Unit tests for the CFC cap module.
671logical function cfc_cap_unit_tests(verbose)
672 logical, intent(in) :: verbose !< If true, output additional
673 !! information for debugging unit tests
674
675 ! Local variables
676 real :: dummy1, dummy2 ! Test values of Schmidt numbers [nondim] or solubilities [mol kg-1 atm-1] for CFC11 and CFC12
677 real :: ta ! A test value of temperature [hectoKelvin]
678 real :: sal ! A test value of salinity [ppt]
679 character(len=120) :: test_name ! Title of the unit test
680
681 cfc_cap_unit_tests = .false.
682 write(stdout,*) '==== MOM_CFC_cap ======================='
683
684 ! test comp_CFC_schmidt, Table 1 in Wanninkhof (2014); doi:10.4319/lom.2014.12.351
685 test_name = 'Schmidt number calculation'
686 call comp_cfc_schmidt(20.0, dummy1, dummy2)
688 compare_values(verbose, test_name, dummy1, 1179.0, 0.5)
690 compare_values(verbose, test_name, dummy2, 1188.0, 0.5)
691
692 if (.not. cfc_cap_unit_tests) write(stdout,'(2x,a)') "Passed "//test_name
693
694 test_name = 'Solubility function, SST = 1.0 C, and SSS = 10 psu'
695 ta = max(0.01, (1.0 + 273.15) * 0.01) ; sal = 10.
696 ! cfc1 = 3.238 10-2 mol kg-1 atm-1
697 ! cfc2 = 7.943 10-3 mol kg-1 atm-1
698 call get_solubility(dummy1, dummy2, ta, sal , 1.0)
700 compare_values(verbose, test_name, dummy1, 3.238e-2, 5.0e-6)
702 compare_values(verbose, test_name, dummy2, 7.943e-3, 5.0e-6)
703
704 if (.not. cfc_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name
705
706 test_name = 'Solubility function, SST = 20.0 C, and SSS = 35 psu'
707 ta = max(0.01, (20.0 + 273.15) * 0.01) ; sal = 35.
708 ! cfc1 = 0.881 10-2 mol kg-1 atm-1
709 ! cfc2 = 2.446 10-3 mol kg-1 atm-1
710 call get_solubility(dummy1, dummy2, ta, sal , 1.0)
712 compare_values(verbose, test_name, dummy1, 8.8145e-3, 5.0e-8)
714 compare_values(verbose, test_name, dummy2, 2.4462e-3, 5.0e-8)
715 if (.not. cfc_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name
716
717end function cfc_cap_unit_tests
718
719!> Test that ans and calc are approximately equal by computing the difference
720!! and comparing it against limit.
721logical function compare_values(verbose, test_name, calc, ans, limit)
722 logical, intent(in) :: verbose !< If true, write results to stdout
723 character(len=80), intent(in) :: test_name !< Brief description of the unit test
724 real, intent(in) :: calc !< computed value in arbitrary units [A]
725 real, intent(in) :: ans !< correct value [A]
726 real, intent(in) :: limit !< value above which test fails [A]
727
728 ! Local variables
729 real :: diff ! Difference in values [A]
730
731 diff = ans - calc
732
733 compare_values = .false.
734 if (diff > limit ) then
735 compare_values = .true.
736 write(stdout,*) "CFC_cap_unit_tests, UNIT TEST FAILED: ", test_name
737 write(stdout,10) calc, ans
738 elseif (verbose) then
739 write(stdout,10) calc, ans
740 endif
741
74210 format("calc=",f22.16," ans",f22.16)
743end function compare_values
744
745!> \namespace mom_CFC_cap
746!!
747!! This module contains the code that is needed to simulate
748!! CFC-11 and CFC-12 using atmospheric and sea ice variables
749!! provided via cap (only NUOPC cap is implemented so far).
750
751end module mom_cfc_cap