source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcoxy_medusa.F90 @ 5707

Last change on this file since 5707 was 5707, checked in by acc, 5 years ago

JPALM —25-08-2015 — add MEDUSA in the branch. MEDUSA version already up-to-date with this trunk revision

File size: 12.0 KB
Line 
1MODULE trcoxy_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcoxy_medusa  ***
4   !! TOP :   MEDUSA
5   !!======================================================================
6   !! History :
7   !!  -   !  2011-07  (A. Yool)             added for ROAM project
8   !!----------------------------------------------------------------------
9#if defined key_medusa && defined key_roam
10   !!----------------------------------------------------------------------
11   !!                                        MEDUSA oxygen cycle
12   !!----------------------------------------------------------------------
13   !!   trc_oxy_medusa        : 
14   !!----------------------------------------------------------------------
15      USE oce_trc
16      USE trc
17      USE sms_medusa
18      USE lbclnk
19      USE prtctl_trc      ! Print control for debugging
20
21      IMPLICIT NONE
22      PRIVATE
23     
24      PUBLIC   trc_oxy_medusa    ! called in trc_bio_medusa
25
26   !!* Substitution
27#  include "domzgr_substitute.h90"
28   !!----------------------------------------------------------------------
29   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
30   !! $Id: trcbio.F90 1146 2008-06-25 11:42:56Z rblod $
31   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
32   !!----------------------------------------------------------------------
33
34! The following is a map of the subroutines contained within this module
35! - trc_oxy_medusa
36!      - CALLS gas_transfer
37!      - CALLS oxy_schmidt
38!      - CALLS oxy_sato
39
40CONTAINS
41
42!=======================================================================
43!
44   SUBROUTINE trc_oxy_medusa( pt, ps, uwind, vwind, pp0, o2, dz, &  !! inputs
45      kw660, o2flux, o2sat )                                        !! outputs
46!     
47!=======================================================================
48      !!
49      !! Title  : Calculates O2 change due to air-sea gas exchange
50      !! Author : Andrew Yool
51      !! Date   : 15/10/04 (revised 08/07/11)
52      !!
53      !! This subroutine calculates oxygen air-sea gas exchange in the
54      !! surface layer of the ocean.  The formulation is taken from one
55      !! written by Ray Najjar for the OCMIP-2 project.  The routine
56      !! calls two other subroutines, oxy_schmidt.f (calculates Schmidt
57      !! number of oxygen) and oxy_sato.f (calculates oxygen saturation
58      !! concentration at 1 atm).
59      !!
60      !! Function inputs are (in order) :
61      !!     pt      temperature                     (degrees C)
62      !!     ps      salinity                        (PSU)
63      !!     uwind   u-wind velocity                 (m/s)
64      !!     vwind   v-wind velocity                 (m/s)
65      !!     pp0     surface pressure                (divided by 1 atm)
66      !!     o2      surface O2 concentration        (mol/m3)
67      !!     dz      surface layer thickness         (m)
68      !! (*) kw660   gas transfer velocity           (m/s)
69      !! (*) o2flux  exchange rate of oxygen         (mol/m3/s)
70      !! (+) o2sat   oxygen saturation concentration (mol/m3)
71      !!
72      !! Where (*) is the function output (note its units). 
73      !!
74!=======================================================================
75
76      implicit none
77!
78      REAL(wp), INTENT( in )    :: pt
79      REAL(wp), INTENT( in )    :: ps
80      REAL(wp), INTENT( in )    :: uwind
81      REAL(wp), INTENT( in )    :: vwind
82      REAL(wp), INTENT( in )    :: pp0
83      REAL(wp), INTENT( in )    :: o2
84      REAL(wp), INTENT( in )    :: dz
85      REAL(wp), INTENT( inout ) :: kw660, o2flux, o2sat
86!
87      REAL(wp) :: scl_wind, kwo2, o2schmidt, o2sato
88!
89! Calculate gas transfer
90!
91      call gas_transfer(uwind, vwind, scl_wind, kw660)
92!
93! Calculate oxygen Schmidt number
94!
95      call oxy_schmidt(pt, o2schmidt)
96!
97! Calculate the transfer velocity for O2 (m/s)
98!
99      kwo2 = kw660 * (660 / o2schmidt)**0.5
100!
101! Calculate the saturation concentration for oxygen (mol/m3)
102!
103      call oxy_sato(pt, ps, o2sato)
104      o2sat = o2sato * pp0
105!
106! Calculate time rate of change of O2 due to gas exchange (mol/m3/s)
107!
108      o2flux = kwo2 * (o2sat - o2) / dz
109!
110      END SUBROUTINE trc_oxy_medusa
111
112!=======================================================================
113!=======================================================================
114!=======================================================================
115
116!=======================================================================
117!
118   SUBROUTINE oxy_schmidt( pt, &  !! input
119      o2_schmidt )                !! output
120!     
121!=======================================================================
122      !!
123      !! Title  : Calculates Schmidt number for ocean uptake of O2
124      !! Author : Andrew Yool
125      !! Date   : 14/10/04 (revised 08/07/11)
126      !!
127      !! This subroutine calculates the Schmidt number for O2 using sea
128      !! surface temperature.  The code is based upon that developed as
129      !! part of the OCMIP-2 project (1998-2000).  The coefficients used
130      !! are taken from Keeling et al. (1998, GBC, 12, 141-163).
131      !!
132      !! Function inputs are (in order) :
133      !!     t           temperature (degrees C)
134      !! (*) o2_schmidt  oxygen Schmidt number
135      !!
136      !! Where (*) is the function output.
137      !!
138!=======================================================================
139
140      implicit none
141!
142      REAL(wp) :: pt, o2_schmidt
143      REAL(wp) :: a0, a1, a2, a3
144!
145      data a0 /    1638.0 /
146      data a1 /    -81.83 /
147      data a2 /     1.483 /
148      data a3 / -0.008004 /
149!
150      o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*a3))
151!
152      END SUBROUTINE oxy_schmidt
153
154!=======================================================================
155!=======================================================================
156!=======================================================================
157
158!=======================================================================
159!
160   SUBROUTINE oxy_sato( pt, ps, &  !! inputs
161      o2_sato )                    !! output
162!     
163!=======================================================================
164      !!
165      !! Title  : Calculates O2 saturation at 1 atm pressure
166      !! Author : Andrew Yool
167      !! Date   : 14/10/04 (revised 08/07/11)
168      !!
169      !! This subroutine calculates the oxygen saturation concentration
170      !! at 1 atmosphere pressure in mol/m3 given ambient temperature
171      !! and salinity.  This formulation is (ostensibly) taken from
172      !! Garcia & Gordon (1992, L&O, 1307-1312).  The function works
173      !! in the range -1.9 <= T <= 40, 0 <= S <= 42.
174      !!
175      !! Function inputs are (in order) :
176      !!     pt       temperature (degrees C)
177      !!     ps       salinity (PSU)
178      !! (*) o2_sato  oxygen saturation (mol/m3)
179      !!
180      !! Where (*) is the function output (note its units). 
181      !!
182      !! Check value : T = 10, S = 35, oxy_sato = 0.282015 mol/m3
183      !!
184!=======================================================================
185
186      implicit none
187!
188      REAL(wp) :: pt, ps, o2_sato
189!
190      REAL(wp) :: a0, a1, a2, a3, a4, a5
191      REAL(wp) :: b0, b1, b2, b3
192      REAL(wp) :: c0
193!
194      REAL(wp) :: tt, tk, ts, ts2, ts3, ts4, ts5
195      REAL(wp) :: ans1, ans2
196!
197      data a0 /  2.00907    /
198      data a1 /  3.22014    /
199      data a2 /  4.05010    /
200      data a3 /  4.94457    /
201      data a4 / -2.56847E-1 /
202      data a5 /  3.88767    /
203!
204      data b0 / -6.24523E-3 /
205      data b1 / -7.37614E-3 /
206      data b2 / -1.03410E-2 /
207      data b3 / -8.17083E-3 /
208!
209      data c0 / -4.88682E-7 /
210!     
211      tt   = 298.15 - pt
212      tk   = 273.15 + pt
213      ts   = log(tt / tk)
214      ts2  = ts**2
215      ts3  = ts**3
216      ts4  = ts**4
217      ts5  = ts**5
218      ans1 = a0 + a1*ts + a2*ts2 + a3*ts3 + a4*ts4 + a5*ts5  &
219           + ps*(b0 + b1*ts + b2*ts2 + b3*ts3)               &
220           + c0*(ps*ps)
221      ans2 = exp(ans1)
222!
223!  Convert from ml/l to mol/m3
224!
225      o2_sato = (ans2 / 22391.6) * 1000.0
226!
227      END SUBROUTINE oxy_sato
228
229!=======================================================================
230!=======================================================================
231!=======================================================================
232
233!=======================================================================
234!
235   SUBROUTINE gas_transfer( uwind, vwind, &  !! input
236      scl_wind, k )                          !! output
237!     
238!=======================================================================
239      !!
240      !! Title  : Calculates gas transfer velocity
241      !! Author : Andrew Yool
242      !! Date   : 15/10/04 (revised 04/08/2011)
243      !!
244      !! This subroutine uses near-surface wind speed to calculate gas
245      !! transfer velocity for use in CO2 and O2 exchange calculations.
246      !!
247      !! Note that the parameterisation of Wanninkhof quoted here is a
248      !! truncation of the original equation.  It excludes a chemical
249      !! enhancement function (based on temperature), although such
250      !! temperature dependence is reported negligible by Etcheto &
251      !! Merlivat (1988).
252      !!
253      !! Note also that in calculating scalar wind, the variance of the
254      !! wind over the period of a timestep is ignored.  Some authors,
255      !! for instance OCMIP-2, favour including some reference to the
256      !! variability of wind.  However, their wind fields are averaged
257      !! over relatively long time periods, and so this issue may be
258      !! safely (!) ignored here.
259      !!
260      !! Subroutine inputs are (in order) :
261      !!     uwind     wind u velocity at 10 m (m/s)
262      !!     vwind     wind v velocity at 10 m (m/s)
263      !! (+) scl_wind  scalar wind velocity at 10 m (m/s)
264      !! (*) k         gas transfer velocity (m/s)
265      !! Where (*) is the function output and (+) is a diagnostic output.
266      !!
267!=======================================================================
268
269      implicit none
270!
271! Input variables
272      REAL(wp) :: uwind, vwind
273!
274! Output variables
275      REAL(wp) :: scl_wind, k, tmp_k
276!
277! Choice of parameterisation
278      INTEGER  :: eqn
279!
280! Coefficients for various parameterisations
281      REAL(wp), DIMENSION(6) :: a
282      REAL(wp), DIMENSION(6) :: b
283!
284! Values of coefficients
285      data a(1) / 0.166 /  ! Liss & Merlivat (1986)    [approximated]
286      data a(2) / 0.3   /  ! Wanninkhof (1992)         [sans enhancement]
287      data a(3) / 0.23  /  ! Nightingale et al. (2000) [good]
288      data a(4) / 0.23  /  ! Nightingale et al. (2000) [better]
289      data a(5) / 0.222 /  ! Nightingale et al. (2000) [best]
290      data a(6) / 0.337 /  ! OCMIP-2                   [sans variability]
291!
292      data b(1) / 0.133 /
293      data b(2) / 0.0   /
294      data b(3) / 0.0   /
295      data b(4) / 0.1   /
296      data b(5) / 0.333 /
297      data b(6) / 0.0   /
298!
299! Which parameterisation is to be used?
300      eqn = 2
301!
302! Calculate scalar wind (m/s)
303      scl_wind = (uwind**2 + vwind**2)**0.5
304!
305! Calculate gas transfer velocity (cm/h)
306      tmp_k = (a(eqn) * scl_wind**2) + (b(eqn) * scl_wind)
307!
308! Convert tmp_k from cm/h to m/s
309      k = tmp_k / (100. * 3600.)
310!
311      END SUBROUTINE gas_transfer
312
313!=======================================================================
314!=======================================================================
315!=======================================================================
316
317#else
318   !!======================================================================
319   !!  Dummy module :                                   No MEDUSA bio-model
320   !!======================================================================
321
322CONTAINS
323
324   SUBROUTINE trc_oxy_medusa( pt, ps, uwind, vwind, pp0, o2, dz, &  !! inputs
325      kw660, o2flux, o2sat )                                        !! outputs
326      USE par_kind
327
328      REAL(wp), INTENT( in )    :: pt
329      REAL(wp), INTENT( in )    :: ps
330      REAL(wp), INTENT( in )    :: pp0
331      REAL(wp), INTENT( in )    :: o2
332      REAL(wp), INTENT( in )    :: dz
333      REAL(wp), INTENT( inout ) :: kw660, o2flux, o2sat
334
335      WRITE(*,*) 'trc_oxy_medusa: You should not have seen this print! error?', kt
336
337   END SUBROUTINE trc_oxy_medusa
338#endif
339
340   !!======================================================================
341END MODULE  trcoxy_medusa
342 
Note: See TracBrowser for help on using the repository browser.