source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcoxy_medusa.F90 @ 9309

Last change on this file since 9309 was 7254, checked in by jpalmier, 4 years ago

JPALM — 17-11-2016 — CMIP6 diags for MEDUSA bugfix

File size: 9.6 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      PUBLIC   oxy_sato          ! called in trc_bio_medusa
26
27   !!* Substitution
28#  include "domzgr_substitute.h90"
29   !!----------------------------------------------------------------------
30   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
31   !! $Id$
32   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34
35! The following is a map of the subroutines contained within this module
36! - trc_oxy_medusa
37!      - CALLS oxy_schmidt
38!      - CALLS oxy_sato
39
40CONTAINS
41
42!=======================================================================
43!
44   SUBROUTINE trc_oxy_medusa( pt, ps, kw660, pp0, o2,  &  !! inputs
45      kwo2, 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      !! AXY (23/06/15): revised to allow common gas transfer velocity
61      !!                 to be used for CO2 and O2; outputs of this
62      !!                 routine amended to mmol/m3 from mol/m3
63      !!
64      !! Function inputs are (in order) :
65      !!     pt      temperature                     (degrees C)
66      !!     ps      salinity                        (PSU)
67      !!     kw660   gas transfer velocity           (m/s)
68      !!     pp0     surface pressure                (divided by 1 atm)
69      !!     o2      surface O2 concentration        (mmol/m3)
70      !! (+) kwo2    gas transfer velocity for O2    (m/s)
71      !! (*) o2flux  exchange rate of oxygen         (mmol/m2/s)
72      !! (+) o2sat   oxygen saturation concentration (mmol/m3)
73      !!
74      !! Where (*) is the function output (note its units). 
75      !!
76!=======================================================================
77
78      implicit none
79!
80      REAL(wp), INTENT( in )    :: pt
81      REAL(wp), INTENT( in )    :: ps
82      REAL(wp), INTENT( in )    :: kw660
83      REAL(wp), INTENT( in )    :: pp0
84      REAL(wp), INTENT( in )    :: o2
85      REAL(wp), INTENT( out )   :: kwo2, o2flux, o2sat
86!
87      REAL(wp) :: o2schmidt, o2sato, mol_o2
88!
89! Oxygen to mol / m3
90!
91      mol_o2 = o2 / 1000.
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 - mol_o2)
109!
110! Oxygen flux and saturation to mmol / m3
111!
112      o2sat  =  o2sat * 1000.
113      o2flux = o2flux * 1000.
114!
115      END SUBROUTINE trc_oxy_medusa
116
117!=======================================================================
118!=======================================================================
119!=======================================================================
120
121!=======================================================================
122!
123   SUBROUTINE oxy_schmidt( pt, &  !! input
124      o2_schmidt )                !! output
125!     
126!=======================================================================
127      !!
128      !! Title  : Calculates Schmidt number for ocean uptake of O2
129      !! Author : Andrew Yool
130      !! Date   : 14/10/04 (revised 08/07/11)
131      !!
132      !! This subroutine calculates the Schmidt number for O2 using sea
133      !! surface temperature.  The code is based upon that developed as
134      !! part of the OCMIP-2 project (1998-2000).  The coefficients used
135      !! are taken from Keeling et al. (1998, GBC, 12, 141-163).
136      !!
137      !! AXY (23/06/2015)
138      !! UPDATED: revised formulation from Wanninkhof (2014) for
139      !! consistency with MOCSY
140      !!
141      !! Winninkhof, R. (2014). Relationship between wind speed and gas
142      !! exchange over the ocean revisited. LIMNOLOGY AND OCEANOGRAPHY-METHODS
143      !! 12, 351-362, doi:10.4319/lom.2014.12.351
144      !!
145      !! Function inputs are (in order) :
146      !!     t           temperature (degrees C)
147      !! (*) o2_schmidt  oxygen Schmidt number
148      !!
149      !! Where (*) is the function output.
150      !!
151!=======================================================================
152
153      implicit none
154!
155      REAL(wp) :: pt, o2_schmidt
156      REAL(wp) :: a0, a1, a2, a3, a4
157!
158! AXY (23/06/15): OCMIP-2 coefficients
159!     data a0 /    1638.0 /
160!     data a1 /    -81.83 /
161!     data a2 /     1.483 /
162!     data a3 / -0.008004 /
163!
164! AXY (23/06/15): Wanninkhof (2014) coefficients
165      data a0 /     1920.4 /
166      data a1 /     -135.6 /
167      data a2 /     5.2121 /
168      data a3 /   -0.10939 /
169      data a4 / 0.00093777 /
170!
171!     o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*a3))
172      o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*(a3 + pt*a4)))
173!
174      END SUBROUTINE oxy_schmidt
175
176!=======================================================================
177!=======================================================================
178!=======================================================================
179
180!=======================================================================
181!
182   SUBROUTINE oxy_sato( pt, ps, &  !! inputs
183      o2_sato )                    !! output
184!     
185!=======================================================================
186      !!
187      !! Title  : Calculates O2 saturation at 1 atm pressure
188      !! Author : Andrew Yool
189      !! Date   : 14/10/04 (revised 08/07/11)
190      !!
191      !! This subroutine calculates the oxygen saturation concentration
192      !! at 1 atmosphere pressure in mol/m3 given ambient temperature
193      !! and salinity.  This formulation is (ostensibly) taken from
194      !! Garcia & Gordon (1992, L&O, 1307-1312).  The function works
195      !! in the range -1.9 <= T <= 40, 0 <= S <= 42.
196      !!
197      !! Function inputs are (in order) :
198      !!     pt       temperature (degrees C)
199      !!     ps       salinity (PSU)
200      !! (*) o2_sato  oxygen saturation (mol/m3)
201      !!
202      !! Where (*) is the function output (note its units). 
203      !!
204      !! Check value : T = 10, S = 35, oxy_sato = 0.282015 mol/m3
205      !!
206!=======================================================================
207
208      implicit none
209!
210      REAL(wp) :: pt, ps, o2_sato
211!
212      REAL(wp) :: a0, a1, a2, a3, a4, a5
213      REAL(wp) :: b0, b1, b2, b3
214      REAL(wp) :: c0
215!
216      REAL(wp) :: tt, tk, ts, ts2, ts3, ts4, ts5
217      REAL(wp) :: ans1, ans2
218!
219      data a0 /  2.00907    /
220      data a1 /  3.22014    /
221      data a2 /  4.05010    /
222      data a3 /  4.94457    /
223      data a4 / -2.56847E-1 /
224      data a5 /  3.88767    /
225!
226      data b0 / -6.24523E-3 /
227      data b1 / -7.37614E-3 /
228      data b2 / -1.03410E-2 /
229      data b3 / -8.17083E-3 /
230!
231      data c0 / -4.88682E-7 /
232!     
233      tt   = 298.15 - pt
234      tk   = 273.15 + pt
235      ts   = log(tt / tk)
236      ts2  = ts**2
237      ts3  = ts**3
238      ts4  = ts**4
239      ts5  = ts**5
240      ans1 = a0 + a1*ts + a2*ts2 + a3*ts3 + a4*ts4 + a5*ts5  &
241           + ps*(b0 + b1*ts + b2*ts2 + b3*ts3)               &
242           + c0*(ps*ps)
243      ans2 = exp(ans1)
244!
245!  Convert from ml/l to mol/m3
246!
247      o2_sato = (ans2 / 22391.6) * 1000.0
248!
249      END SUBROUTINE oxy_sato
250
251!=======================================================================
252!=======================================================================
253!=======================================================================
254
255#else
256   !!======================================================================
257   !!  Dummy module :                                   No MEDUSA bio-model
258   !!======================================================================
259
260CONTAINS
261
262   SUBROUTINE trc_oxy_medusa( pt, ps, kw660, pp0, o2,  &  !! inputs
263      o2flux, o2sat )                                     !! outputs
264      USE par_kind
265
266      REAL(wp), INTENT( in )    :: pt
267      REAL(wp), INTENT( in )    :: ps
268      REAL(wp), INTENT( in )    :: kw660
269      REAL(wp), INTENT( in )    :: pp0
270      REAL(wp), INTENT( in )    :: o2
271      REAL(wp), INTENT( inout ) :: o2flux, o2sat
272
273      WRITE(*,*) 'trc_oxy_medusa: You should not have seen this print! error?', kt
274
275   END SUBROUTINE trc_oxy_medusa
276#endif
277
278   !!======================================================================
279END MODULE  trcoxy_medusa
280 
Note: See TracBrowser for help on using the repository browser.