New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
mocsy_wrapper.F90 in branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_wrapper.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.

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