forked from AquaticEcoDynamics/libaed2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
aed2_oxygen.F90
410 lines (342 loc) · 19.4 KB
/
aed2_oxygen.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
!###############################################################################
! !
! .----------------. .----------------. .----------------. !
! | .--------------. || .--------------. || .--------------. | !
! | | ____ | || | ____ ____ | || | ____ ____ | | !
! | | .' `. | || | |_ _||_ _| | || | |_ _||_ _| | | !
! | | / .--. \ | || | \ \ / / | || | \ \ / / | | !
! | | | | | | | || | > `' < | || | \ \/ / | | !
! | | \ `--' / | || | _/ /'`\ \_ | || | _| |_ | | !
! | | `.____.' | || | |____||____| | || | |______| | | !
! | | | || | | || | | | !
! | '--------------' || '--------------' || '--------------' | !
! '----------------' '----------------' '----------------' !
! !
!###############################################################################
!# #
!# aed2_oxygen.F90 #
!# #
!# Developed by : #
!# AquaticEcoDynamics (AED) Group #
!# The University of Western Australia #
!# #
!# http://aquatic.science.uwa.edu.au/ #
!# #
!# Copyright 2013 - 2019 - The University of Western Australia #
!# #
!# GLM is free software: you can redistribute it and/or modify #
!# it under the terms of the GNU General Public License as published by #
!# the Free Software Foundation, either version 3 of the License, or #
!# (at your option) any later version. #
!# #
!# GLM is distributed in the hope that it will be useful, #
!# but WITHOUT ANY WARRANTY; without even the implied warranty of #
!# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
!# GNU General Public License for more details. #
!# #
!# You should have received a copy of the GNU General Public License #
!# along with this program. If not, see <http://www.gnu.org/licenses/>. #
!# #
!# ----------------------------------------------------------------------- #
!# #
!# Created May 2011 #
!# Track changes on GitHub @ https://github.com/AquaticEcoDynamics/libaed2 #
!# #
!###############################################################################
#include "aed2.h"
MODULE aed2_oxygen
!-------------------------------------------------------------------------------
! aed2_oxygen --- oxygen biogeochemical model
!
! The AED module oxygen contains equations that describe exchange of
! oxygen across the air/water interface and sediment flux. Other modules can
! also add or consume oxygen. A summary fo the module is provided online at:
! http://aquatic.science.uwa.edu.au/research/models/AED/aed_oxygen.html
!
!-------------------------------------------------------------------------------
USE aed2_core
USE aed2_util, ONLY: aed2_gas_piston_velocity, aed2_oxygen_sat
IMPLICIT NONE
PRIVATE
!
PUBLIC aed2_oxygen_data_t
!
TYPE,extends(aed2_model_data_t) :: aed2_oxygen_data_t
!# Variable identifiers
INTEGER :: id_oxy
INTEGER :: id_Fsed_oxy
INTEGER :: id_oxy_sat
INTEGER :: id_atm_oxy_exch, id_sed_oxy
INTEGER :: id_sed_oxy_pel, id_atm_oxy_exch3d
INTEGER :: id_temp, id_salt, id_wind
INTEGER :: id_larea, id_lht, id_cell_vel
!# Model parameters
AED_REAL :: Fsed_oxy,Ksed_oxy,theta_sed_oxy
INTEGER :: oxy_piston_model
LOGICAL :: use_sed_model
CONTAINS
PROCEDURE :: define => aed2_define_oxygen
PROCEDURE :: calculate_surface => aed2_calculate_surface_oxygen
PROCEDURE :: calculate => aed2_calculate_oxygen
PROCEDURE :: calculate_benthic => aed2_calculate_benthic_oxygen
! PROCEDURE :: mobility => aed2_mobility_oxygen
! PROCEDURE :: light_extinction => aed2_light_extinction_oxygen
! PROCEDURE :: particle_bgc => aed2_particle_bgc_oxygen
! PROCEDURE :: delete => aed2_delete_oxygen
END TYPE
! MODULE GLOBALS
INTEGER :: diag_level = 10 ! 0 = no diagnostic outputs
! 1 = basic diagnostic outputs
! 2-10 = most diagnostic outputs
! >10 = debug/checking outputs
!===============================================================================
CONTAINS
!###############################################################################
SUBROUTINE aed2_define_oxygen(data, namlst)
!-------------------------------------------------------------------------------
! Setup and initialise the aed2_oxygen model
!
! Here, the oxygen namelist is read and te variables exported
! by the model are registered with AED2.
!-------------------------------------------------------------------------------
!ARGUMENTS
CLASS (aed2_oxygen_data_t),INTENT(inout) :: data
INTEGER,INTENT(in) :: namlst
!
!LOCALS
INTEGER :: status
! %% NAMELIST
AED_REAL :: oxy_initial = 300. !% initial dissolved oxygen (DO) concentration
!% $$mmol\,m^{-3}$$
!% float
!%
!% 0 - 400
!% Note: will be overwritten by GLM or TFV IC
AED_REAL :: oxy_min = 0. !% minimum dissolved oxygen (DO) concentration
!% $$mmol\,m^{-3}$$
!% float
!%
!%
!% Optional variable to enforce negative number clipping
AED_REAL :: oxy_max = nan_ !% maxmium dissolved oxygen (DO) concentration
!% $$mmol\,m^{-3}$$
!% float
!% -
!% 1000
!% Optional variable to enforce high number clipping
AED_REAL :: Fsed_oxy = -20.0 !% sediment oxygen demand (SOD)
!% $$mmol\,m^{-2}\,day^{-1}$$
!% float
!%
!% -100
!% Note: unused if Fsed_oxy_variable is activated via aed2_sedflux
AED_REAL :: Ksed_oxy = 30.0 !% half-saturation concentration of oxygen sediment flux
!% $$mmol\,m^{-3}$$
!% float
!% 50
!% 10-100
!% Changes the sensitivity of the oxygen flux to the
!- overlying oxygen concentration
AED_REAL :: theta_sed_oxy = 1.0 !% Arrhenius temperature multiplier for sediment oxygen flux
!% -
!% float
!% 1e+00
!% 1.04 - 1.12
!% Changes the sensitivity of the oxygen flux to
!- the overlying temperature
CHARACTER(len=64) :: Fsed_oxy_variable ='' !% oxygen sediment flux variable link
!% -
!% string
!% '
!% e.g.: SDF_Fsed_oxy
!% will use the value supplied by the aed2_sedflux
!- model for Fsed_oxy; use this option to allow for
!- spatial or temperal variation
INTEGER :: oxy_piston_model = 1 !% specifies the atm exchange piston model
!% -
!% integer
!% 1
!% 1 - X
!% Choice depends on waterbody type
! %% END NAMELIST
NAMELIST /aed2_oxygen/ oxy_initial, oxy_min, oxy_max, &
Fsed_oxy, Ksed_oxy, theta_sed_oxy, &
Fsed_oxy_variable, oxy_piston_model, &
diag_level
!
!-------------------------------------------------------------------------------
!BEGIN
print *," aed2_oxygen initialization"
! Read the namelist
read(namlst,nml=aed2_oxygen,iostat=status)
IF (status /= 0) STOP 'Error reading namelist aed2_oxygen'
! Store parameter values in the modules own derived type
! NB: all rates must be provided in values per day,
! and are converted here to values per second.
data%Ksed_oxy = Ksed_oxy
data%Fsed_oxy = Fsed_oxy/secs_per_day
data%theta_sed_oxy = theta_sed_oxy
data%use_sed_model = Fsed_oxy_variable .NE. ''
data%oxy_piston_model = oxy_piston_model
! Register state variables
data%id_oxy = aed2_define_variable('oxy','mmol/m**3','oxygen', &
oxy_initial,minimum=oxy_min,maximum=oxy_max)
! Register the link to external variables
IF (data%use_sed_model) data%id_Fsed_oxy = aed2_locate_global_sheet(Fsed_oxy_variable)
! Register diagnostic variables
IF (diag_level>0) THEN
data%id_oxy_sat = aed2_define_diag_variable( &
'sat', '%', 'oxygen saturation')
data%id_sed_oxy = aed2_define_sheet_diag_variable( &
'sed_oxy', 'mmol/m**2/d', 'O2 exchange across sed/water interface')
data%id_atm_oxy_exch = aed2_define_sheet_diag_variable( &
'atm_oxy_flux', 'mmol/m**2/d', 'O2 exchange across atm/water interface')
IF (diag_level>10) THEN
data%id_sed_oxy_pel = aed2_define_diag_variable( &
'sed_oxy_pel', 'mmol/m**2/d', 'O2 exchange across sed/water interface')
data%id_atm_oxy_exch3d = aed2_define_diag_variable( &
'atm_oxy_exch3d', 'mmol/m**3/d', 'Oxygen exchange across atm/water interface')
ENDIF
ENDIF
! Register environmental dependencies
data%id_temp = aed2_locate_global('temperature') ! Temperature (degrees Celsius)
data%id_salt = aed2_locate_global('salinity') ! Salinity (psu)
! data%id_pres = aed2_locate_global_sheet('pressure') ! Pressure (dbar = 10 kPa)
data%id_wind = aed2_locate_global_sheet('wind_speed') ! Wind speed at 10 m above surface (m/s)
data%id_larea = aed2_locate_global_sheet('layer_area')
data%id_lht = aed2_locate_global('layer_ht')
data%id_cell_vel = -1
IF( oxy_piston_model>3 )data%id_cell_vel= aed2_locate_global('cell_vel')! needed for k600
END SUBROUTINE aed2_define_oxygen
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!###############################################################################
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!###############################################################################
SUBROUTINE aed2_calculate_surface_oxygen(data,column,layer_idx)
!-------------------------------------------------------------------------------
! Air-water exchange for the aed oxygen model
!-------------------------------------------------------------------------------
!ARGUMENTS
CLASS (aed2_oxygen_data_t),INTENT(in) :: data
TYPE (aed2_column_t),INTENT(inout) :: column(:)
INTEGER,INTENT(in) :: layer_idx
!
!LOCALS
! Environment
AED_REAL :: temp, salt, wind, depth
AED_REAL :: vel = 0.0001
! State
AED_REAL :: oxy
! Temporary variables
AED_REAL :: oxy_atm_flux = zero_ ! Surface atm flux of O2
AED_REAL :: Coxy_air = zero_ ! Dissolved oxygen in the air phase
AED_REAL :: koxy_trans = zero_ ! k600 for O2
AED_REAL :: windHt = 10.0 ! Height of U10 sensor
AED_REAL :: f_pres = 1.0 ! Pressure correction, only applicable at high altitudes
!
!-------------------------------------------------------------------------------
!BEGIN
!Get dependent state variables from physical driver
temp = _STATE_VAR_(data%id_temp) ! Temperature (degrees Celsius)
salt = _STATE_VAR_(data%id_salt) ! Salinity (psu)
wind = _STATE_VAR_S_(data%id_wind) ! Wind speed at 10 m above surface (m/s)
windHt= 10. ! Assumed wind height of 10m
depth = MAX( _STATE_VAR_(data%id_lht), one_ )
IF (data%id_cell_vel > 0 ) vel = _STATE_VAR_(data%id_cell_vel)
! Retrieve current (local) state variable values.
oxy = _STATE_VAR_(data%id_oxy) ! Concentration of oxygen in surface layer
!koxy_trans = aed2_gas_piston_velocity(windHt,wind,temp,salt)
koxy_trans = aed2_gas_piston_velocity(windHt,wind,temp,salt, &
vel=vel, &
depth=depth, &
schmidt_model=2, &
piston_model=data%oxy_piston_model)
! First get the oxygen concentration in the air phase at the interface
! (taken from Riley and Skirrow, 1974)
f_pres = 1.0 ! set pressure function here using data%id_pres
Coxy_air = f_pres * aed2_oxygen_sat(salt,temp)
! Get the oxygen flux
oxy_atm_flux = koxy_trans * (Coxy_air - oxy)
! Transfer surface exchange value to AED2 (mmmol/m2/s) converted by driver
_FLUX_VAR_T_(data%id_oxy) = oxy_atm_flux
! Also store oxygen flux across the atm/water interface as a diagnostic (mmmol/m2/day)
IF (diag_level>0) THEN
_DIAG_VAR_S_(data%id_atm_oxy_exch) = oxy_atm_flux * secs_per_day
_DIAG_VAR_(data%id_oxy_sat) = Coxy_air
ENDIF
! Also store oxygen flux across the atm/water interface as a diagnostic (mmmol/m2/day)
IF (diag_level>10) &
_DIAG_VAR_(data%id_atm_oxy_exch3d) = oxy_atm_flux * secs_per_day / _STATE_VAR_(data%id_lht)
END SUBROUTINE aed2_calculate_surface_oxygen
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!###############################################################################
SUBROUTINE aed2_calculate_oxygen(data,column,layer_idx)
!-------------------------------------------------------------------------------
! Right hand sides of aed2_oxygen model
!-------------------------------------------------------------------------------
!ARGUMENTS
CLASS(aed2_oxygen_data_t),INTENT(in) :: data
TYPE (aed2_column_t),INTENT(inout) :: column(:)
INTEGER,INTENT(in) :: layer_idx
!
!LOCALS
AED_REAL :: oxy, temp, salt
AED_REAL :: f_pres, coxy_sat
!-------------------------------------------------------------------------------
!BEGIN
! Get dependent state variables from physical driver
temp = _STATE_VAR_(data%id_temp) ! Temperature (degrees Celsius)
salt = _STATE_VAR_(data%id_salt) ! Salinity (psu)
! Retrieve current (local) state variable values.
oxy = _STATE_VAR_(data%id_oxy)! oxygen
! Compute the oxygen saturation for diagnostic output
f_pres = 1.0 ! set pressure function here using data%id_pres
coxy_sat = f_pres * aed2_oxygen_sat(salt,temp)
! Export diagnostic variables
_DIAG_VAR_(data%id_oxy_sat) = (oxy/coxy_sat)*100.
END SUBROUTINE aed2_calculate_oxygen
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!###############################################################################
SUBROUTINE aed2_calculate_benthic_oxygen(data,column,layer_idx)
!-------------------------------------------------------------------------------
! Calculate pelagic bottom fluxes and benthic sink and source terms of AED oxygen.
! Everything in units per surface area (not volume!) per time.
!-------------------------------------------------------------------------------
!ARGUMENTS
CLASS (aed2_oxygen_data_t),INTENT(in) :: data
TYPE (aed2_column_t),INTENT(inout) :: column(:)
INTEGER,INTENT(in) :: layer_idx
!
!LOCALS
! Environment
AED_REAL :: temp !, layer_ht
! State
AED_REAL :: oxy
! Temporary variables
AED_REAL :: oxy_flux, Fsed_oxy
!
!-------------------------------------------------------------------------------
!BEGIN
! Retrieve current environmental conditions for the bottom pelagic layer.
temp = _STATE_VAR_(data%id_temp) ! local temperature
! Retrieve current (local) state variable values.
oxy = _STATE_VAR_(data%id_oxy)! oxygen
IF (data%use_sed_model) THEN
Fsed_oxy = _DIAG_VAR_S_(data%id_Fsed_oxy)
ELSE
Fsed_oxy = data%Fsed_oxy
ENDIF
! Compute the sediment flux dependent on overlying oxygen & temperature
oxy_flux = Fsed_oxy * MIN(3.,oxy/(data%Ksed_oxy+oxy) * (data%theta_sed_oxy**(temp-20.0)))
! Set bottom fluxes for the pelagic (change per surface area per second)
! Transfer sediment flux value to AED2
_FLUX_VAR_(data%id_oxy) = _FLUX_VAR_(data%id_oxy) + (oxy_flux)
! Set sink and source terms for the benthos (change per surface area per second)
! Note that this must include the fluxes to and from the pelagic.
!_FLUX_VAR_B_(data%id_ben_oxy) = _FLUX_VAR_B_(data%id_ben_oxy) + (-oxy_flux)
! Also store sediment flux as diagnostic variable.
_DIAG_VAR_S_(data%id_sed_oxy) = oxy_flux * secs_per_day
IF (diag_level>10) _DIAG_VAR_(data%id_sed_oxy_pel) = oxy_flux * secs_per_day
END SUBROUTINE aed2_calculate_benthic_oxygen
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
END MODULE aed2_oxygen