source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_wrapper.F90 @ 7611

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

JPALM — 26-01-2017 — Fix MOCSY depth0 init value

File size: 14.3 KB
Line 
1MODULE mocsy_wrapper
2   !!======================================================================
3   !!                         ***  MODULE mocsy_wrapper  ***
4   !! TOP :   MEDUSA
5   !!======================================================================
6   !! History :
7   !!  -   !  2015-06  (A. Yool)             added for UKESM project
8   !!----------------------------------------------------------------------
9#if defined key_medusa && defined key_roam && defined key_mocsy
10   !!----------------------------------------------------------------------
11   !!                                        MEDUSA carbonate chemistry
12   !!----------------------------------------------------------------------
13   !!   trc_co2_medusa        : 
14   !!----------------------------------------------------------------------
15
16   USE mocsy_mainmod
17   USE mocsy_gasflux
18   USE in_out_manager  ! I/O manager
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC mocsy_interface  ! called in trc_bio_medusa
24   PUBLIC mocsy_carbchem   ! called in trc_bio_medusa
25
26   !!* Substitution
27#  include "domzgr_substitute.h90"
28   !!----------------------------------------------------------------------
29   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
30   !! $Id$
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! - mocsy_interface
36!    - CALLS mocsy_carbchem
37!
38! - mocsy_carbchem
39!    - CALLS vars
40!      - CALLS p80
41!      - CALLS constants
42!      - CALLS rho
43!      - CALLS varsolver
44!        - CALLS [phsolvers routines]
45!      - CALLS varsolver (again)
46!    - CALLS x2pCO2atm
47!    - CALLS p2fCO2
48!    - CALLS schmidt_co2
49!    - CALLS surface_K0
50
51CONTAINS
52
53!=======================================================================
54!
55      SUBROUTINE mocsy_interface( temp, sal, alk, dic, sil, phos, &
56           Patm, depth, lat, kw660, xco2, N,                      &
57           ph, pco2, fco2, co2, hco3, co3, OmegaA,                &
58           OmegaC, BetaD, rhoSW, p, tempis,                       &
59           pco2atm, fco2atm, schmidtco2, kwco2, K0,               &
60           co2starair, co2flux, dpco2 )
61!
62!=======================================================================
63!
64! AXY (26/06/15): to preserve both MEDUSA's scalar variables and
65!                 MOCSY's vector variables on code traceability
66!                 grounds, this additional wrapper does a rather
67!                 superfluous conversion between the two; in the
68!                 fullness of time, this should be dispensed with
69!                 and MOCSY used in its full vector form
70!
71        USE mocsy_singledouble
72        IMPLICIT NONE
73
74!> ======================================================================
75!  VARIABLES
76!> ======================================================================
77!
78   INTEGER, INTENT(in) :: N
79!
80!  MEDUSA-side
81!  Input variables
82   REAL(kind=wp), INTENT(in)  :: temp, sal, alk, dic, sil, phos
83   REAL(kind=wp), INTENT(in)  :: Patm, depth, lat, kw660, xco2
84!
85!  Output variables
86   REAL(kind=wp), INTENT(out) :: ph, pco2, fco2, co2, hco3, co3, OmegaA
87   REAL(kind=wp), INTENT(out) :: OmegaC, BetaD, rhoSW, p, tempis
88   REAL(kind=wp), INTENT(out) :: pco2atm, fco2atm, schmidtco2, kwco2, K0
89   REAL(kind=wp), INTENT(out) :: co2starair, co2flux, dpco2
90!
91!  MOCSY-side
92!  Input variables
93   REAL(kind=wp), DIMENSION(N) :: mtemp, msal, malk, mdic, msil, mphos
94   REAL(kind=wp), DIMENSION(N) :: mPatm, mdepth, mlat, mkw660, mxco2
95!
96!  Output variables
97   REAL(kind=wp), DIMENSION(N) :: mph, mpco2, mfco2, mco2, mhco3, mco3, mOmegaA
98   REAL(kind=wp), DIMENSION(N) :: mOmegaC, mBetaD, mrhoSW, mp, mtempis
99   REAL(kind=wp), DIMENSION(N) :: mpco2atm, mfco2atm, mschmidtco2, mkwco2, mK0
100   REAL(kind=wp), DIMENSION(N) :: mco2starair, mco2flux, mdpco2
101!
102!> ----------------------------------------------------------------------
103!  Set MOCSY inputs to equal MEDUSA inputs (amend units here)
104!> ----------------------------------------------------------------------
105!
106      mtemp(1)  = temp         ! degrees C
107      msal(1)   = sal          ! PSU
108      malk(1)   = alk / 1000.  ! meq  / m3 -> eq  / m3
109      mdic(1)   = dic / 1000.  ! mmol / m3 -> mol / m3
110      msil(1)   = sil / 1000.  ! mmol / m3 -> mol / m3
111      mphos(1)  = phos / 1000. ! mmol / m3 -> mol / m3
112      mPatm(1)  = Patm         ! atm
113      mdepth(1) = depth        ! m
114      mlat(1)   = lat          ! degrees N
115      mkw660(1) = kw660        ! m / s
116      mxco2(1)  = xco2         ! ppm
117!
118!> ----------------------------------------------------------------------
119!  Call MOCSY
120!> ----------------------------------------------------------------------
121!
122      CALL mocsy_carbchem( mtemp, msal, malk, mdic, msil, mphos, & ! inputs
123      mPatm, mdepth, mlat, mkw660, mxco2, 1,                     & ! inputs
124      mph, mpco2, mfco2, mco2, mhco3, mco3, mOmegaA,             & ! outputs
125      mOmegaC, mBetaD, mrhoSW, mp, mtempis,                      & ! outputs
126      mpco2atm, mfco2atm, mschmidtco2, mkwco2, mK0,              & ! outputs
127      mco2starair, mco2flux, mdpco2 ) ! outputs
128!
129!> ----------------------------------------------------------------------
130!  Set MOCSY outputs to equal MEDUSA outputs (amend units here)
131!> ----------------------------------------------------------------------
132!
133      ph         = mph(1)                 ! standard units
134      pco2       = mpco2(1)               ! uatm
135      fco2       = mfco2(1)               ! uatm
136      co2        = mco2(1) * 1000.        ! mol / m3 -> mmol / m3
137      hco3       = mhco3(1) * 1000.       ! mol / m3 -> mmol / m3
138      co3        = mco3(1) * 1000.        ! mol / m3 -> mmol / m3
139      OmegaA     = mOmegaA(1)             ! dimensionless
140      OmegaC     = mOmegaC(1)             ! dimensionless
141      BetaD      = mBetaD(1)              ! dimensionless
142      rhoSW      = mrhoSW(1)              ! kg / m3
143      p          = mp(1)                  ! db
144      tempis     = mtempis(1)             ! degrees C
145      pco2atm    = mpco2atm(1)            ! uatm
146      fco2atm    = mfco2atm(1)            ! uatm
147      schmidtco2 = mschmidtco2(1)         ! dimensionless
148      kwco2      = mkwco2(1)              ! m / s
149      K0         = mK0(1)                 ! (mol/kg) / atm
150      co2starair = mco2starair(1) * 1000. ! mol / m3 -> mmol / m3
151      co2flux    = mco2flux(1) * 1000.    ! mol / m2 / s -> mmol / m2 / s
152      dpco2      = mdpco2(1)              ! uatm
153
154  RETURN
155
156  END SUBROUTINE
157
158!-----------------------------------------------------------------------
159
160!=======================================================================
161!
162      SUBROUTINE mocsy_carbchem( temp, sal, alk, dic, sil, phos,  &
163           Patm, depth, lat, kw660, xco2, N,                      &
164           ph, pco2, fco2, co2, hco3, co3, OmegaA,                &
165           OmegaC, BetaD, rhoSW, p, tempis,                       &
166           pco2atm, fco2atm, schmidtco2, kwco2, K0,               &
167           co2starair, co2flux, dpco2 )
168!
169!=======================================================================
170!
171! AXY (23/06/15): MOCSY introduced to MEDUSA to update carbonate
172!                 chemistry for UKESM1 project
173!
174! Orr, J. C. and Epitalon, J.-M.: Improved routines to model the
175! ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485-499,
176! doi:10.5194/gmd-8-485-2015, 2015.
177!
178! "mocsy is a Fortran 95 package designed to compute all ocean
179! carbonate system variables from DIC and total Alk, particularly
180! from models.  It updates previous OCMIP code, avoids 3 common
181! model approximations, and offers the best-practice constants as
182! well as more recent options. Its results agree with those from
183! CO2SYS-MATLAB within 0.005%."
184!
185! Where possible the code remains identical to that published by
186! Orr & Epitalon (2015; henceforth OE15); some consolidation of
187! MOCSY files has taken place, and the resulting package is
188! comprised of four modules:
189!                 
190! 1. mocsy_wrapped.f90
191!    - This module contains the interface between MEDUSA and the
192!      main MOCSY routines; it draws from (but is very different
193!      from) the test_mocsy.f90 routine provided by OE15
194!
195! 2. mocsy_singledouble.f90
196!    - This module contains only precision definitions; it is
197!      based on the singledouble.f90 routine provided by OE15
198!
199! 3. mocsy_phsolvers.f90
200!    - This module contains a suite of solvers derived from
201!      Munhoven (2013) and modified by OE15; it is based on the
202!      phsolver.f90 routine provided by OE15
203!
204! 4. mocsy_gasflux.f90
205!    - This module contains a series of subroutines used to
206!      calculate the air-sea flux of CO2; it consolidates four
207!      MOCSY routines with an OCMIP-2 Schmidt number routine
208!      updated with the new parameterisations of Wanninkhof (2014)
209!
210! 5. mocsy_mainmod.f90
211!    - This module consolidates all of the remaining MOCSY
212!      functions and subroutines from OE15 into a single file
213!      for convenience
214!
215! NOTE: it is still possible to run PML's carbonate chemistry
216! routine in MEDUSA; at present this remains the default, if
217! less-preferred, routine (i.e. is used if key_mocsy is not
218! present)
219!
220        USE mocsy_singledouble
221        IMPLICIT NONE
222
223!> ======================================================================
224!  VARIABLES
225!> ======================================================================
226!
227!  Input variables
228   REAL(kind=wp), INTENT(in),  DIMENSION(N) :: temp, sal, alk, dic, sil, phos
229   REAL(kind=wp), INTENT(in),  DIMENSION(N) :: Patm, depth, lat, kw660, xco2
230   INTEGER, INTENT(in) :: N
231!
232!  Output variables
233   REAL(kind=wp), INTENT(out), DIMENSION(N) :: ph, pco2, fco2, co2, hco3, co3, OmegaA
234   REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaC, BetaD, rhoSW, p, tempis
235   REAL(kind=wp), INTENT(out), DIMENSION(N) :: pco2atm, fco2atm, schmidtco2, kwco2, K0
236   REAL(kind=wp), INTENT(out), DIMENSION(N) :: co2starair, co2flux, dpco2
237!
238!  Local variables
239   REAL(kind=wp), DIMENSION(N) :: depth0, co2star
240   INTEGER :: i, totdat
241!
242!  "vars" Input options
243   CHARACTER(10) :: optCON, optT, optP, optB, optKf, optK1K2
244!   
245!  initialise depth0 to 0
246   depth0 = 0.0
247   
248!> ======================================================================
249!  CONFIGURE OPTIONS
250!> ======================================================================
251!> OPTIONS: see complete documentation in 'vars' subroutine
252!> Typical options for MODELS
253   optCON  = 'mol/m3' ! input concentrations are in MOL/M3
254   optT    = 'Tpot'   ! input temperature, variable 'temp', is POTENTIAL temp [°C]
255   optP    = 'm'      ! input variable 'depth' is in METERS
256   optB    = 'l10'    ! Lee et al. (2010) formulation for total boron
257   optK1K2 = 'l'      ! Lueker et al. (2000) formulations for K1 & K2 (best practices)
258   optKf   = 'dg'     ! Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994)
259
260!> ======================================================================
261!  CARBONATE CHEMISTRY CALCULATIONS
262!> ======================================================================
263!> Call mocsy's main subroutine to compute carbonate system variables:
264!> pH, pCO2, fCO2, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R
265!> FROM temperature, salinity, total alkalinity, dissolved inorganic
266!> carbon, silica, phosphate, depth (or pressure) (1-D arrays)
267   call vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis,  &  ! OUTPUT
268             temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,                      &  ! INPUT
269             optCON, optT, optP)
270!            optCON, optT, optP, optB, optK1K2, optKf)                                    ! INPUT OPTIONS
271
272!> ======================================================================
273!  GAS EXCHANGE CALCULATIONS
274!> ======================================================================
275!>
276!> Only calculate gas exchange fields if depth = 0
277   if (depth(1) .eq. 0.0) then
278!
279!     Compute pCO2atm [uatm] from xCO2 [ppm], atmospheric pressure [atm], & vapor pressure of seawater
280!     pCO2atm = (Patm - pH20(i)) * xCO2,   where pH20 is the vapor pressure of seawater [atm]
281      CALL x2pCO2atm(xco2, temp, sal, Patm, N, & ! INPUT
282      pco2atm)                                   ! OUTPUT
283
284!     Compute fCO2atm [uatm] from pCO2atm [uatm] & fugacity coefficient [unitless]
285!     fCO2atm = pCO2atm * fugcoeff,   where fugcoeff= exp(Patm*(B + 2.0*xc2*Del)/(R*tk) )
286      CALL p2fCO2(pco2atm, temp, Patm, depth0, N, & ! INPUT
287      fco2atm)                                      ! OUTPUT
288
289!     Compute Schmidt number for CO2 from potential temperature
290      CALL schmidt_co2(temp, N, & ! INPUT
291      schmidtco2)                 ! OUTPUT
292
293!     Compute transfer velocity for CO2 in m/s (see equation [4] in OCMIP2 design document & OCMIP2 Abiotic HOWTO)
294      kwco2 = kw660 * (660./schmidtco2)**0.5
295
296!     Surface K0 [(mol/kg) / atm] at T, S of surface water
297      CALL surface_K0(temp, sal, N, & ! INPUT
298      K0)                             ! OUTPUT
299
300!     "Atmospheric" [CO2*], air-sea CO2 flux, sfc DIC rate of change, & Delta pCO2
301!     all "lifted" from the gasx.f90 function of MOCSY
302      co2starair = K0 * fco2atm * 1.0e-6_wp * rhoSW ! Equil. [CO2*] for atm CO2 at Patm & sfc-water T,S [mol/m3]
303      co2star    = co2                              ! Oceanic [CO2*] in [mol/m3] from vars.f90
304      co2flux    = kwco2 * (co2starair - co2star)   ! Air-sea CO2 flux [mol/(m2 * s)]
305!     co2ex      = co2flux / dz1                    ! Change in sfc DIC due to gas exchange [mol/[m3 * s)]
306      dpco2      = pco2 - pco2atm                   ! Delta pCO2 (oceanic - atmospheric pCO2) [uatm]
307!     
308   endif
309
310  RETURN
311
312  END SUBROUTINE
313
314!-----------------------------------------------------------------------
315
316!=======================================================================
317!=======================================================================
318!=======================================================================
319
320#  else
321   !!======================================================================
322   !!  Dummy module :                         No MOCSY carbonate chemistry
323   !!======================================================================
324
325CONTAINS
326
327   SUBROUTINE mocsy_interface( kt )                   ! Empty routine
328
329      INTEGER, INTENT( in ) ::   kt
330
331      WRITE(*,*) 'mocsy_interface: You should not have seen this print! error?', kt
332
333   END SUBROUTINE mocsy_interface
334#endif 
335
336   !!======================================================================
337
338END MODULE mocsy_wrapper
Note: See TracBrowser for help on using the repository browser.