source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_wrapper.F90

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

JPALM — 11-04-2017 — MEDUSA spring tidy-up refreshning session

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