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.
h3cflx.F in trunk/NEMO/TOP_SRC/SMS – NEMO

source: trunk/NEMO/TOP_SRC/SMS/h3cflx.F @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.6 KB
Line 
1
2CCC $Header$ 
3CCC  TOP 1.0 , LOCEAN-IPSL (2005) 
4C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
5C ---------------------------------------------------------------------------
6CCC $Header$
7CDIR$ LIST
8      SUBROUTINE h3cflx
9#if defined key_passivetrc && defined key_trc_hamocc3
10CCC---------------------------------------------------------------------
11CCC
12CCC                       ROUTINE h3cflx
13CCC                     ******************
14CCC
15CCC
16CC     PURPOSE.
17CC     --------
18CC          *H3CFLX* CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
19CC
20CC     METHOD.
21CC     -------
22CC          SOLVING SYSTEM OF TWO NON-LINEAR SIMULTANEOUS EQUATIONS
23CC     FOR [H2CO3] AND [H+] 
24CC
25CC     EXTERNALS.
26CC     ----------
27CC          NONE.
28CC
29CC     REFERENCE.
30CC     ----------
31CC
32CC          BROECKER, W.S., AND T.-H. PENG (1982)
33CC          TRACERS IN THE SEA.
34CC          ELDIGIO PRESS, LAMONT-DOHERTY GEOLOGICAL OBSERVATORY,
35CC          PALISADES, N.Y., 690 PP..
36CC
37CC          MOOK, W.G. (1986)
38CC          13C IN ATMOSPHERIC CO2.
39CC          NETHERLANDS JOURNAL OF SEA RESEARCH, 20(2/3): 211-223.
40CC
41CC          SCARBOROUGH, J. (1958) NUMERICAL MATHEMATICAL ANALYSIS.
42CC          OXFORD UNIVERSITY PRESS, LONDON, 4TH ED., 576 PP..
43CC
44CC*     VARIABLE           TYPE    PURPOSE.
45CC      --------           ----    --------
46CC
47CC      *ZEROAB*           REAL    273.15 (0 DEG C IN DEG K), HALF PREC.
48CC      *ONR*              REAL    1., HALF PRECISION
49CC      *TWR*              REAL    2., HALF PRECISION
50CC      *FOURR*            REAL    4., HALF PRECISION
51CC      *EVFR00*           REAL    COEFFICIENT OF TEMPERATURE DEPENDENT
52CC                                 ISOTOPIC CARBON FRACTIONATION FACTOR
53CC                                 FOR EVAPORATION (MOOK, 1986)
54CC      *EVFR01*           REAL    COEFFICIENT OF TEMPERATURE DEPENDENT
55CC                                 ISOTOPIC CARBON FRACTIONATION FACTOR
56CC                                 FOR EVAPORATION (MOOK, 1986)
57CC      *XN*               REAL    SURFACE [H+] (EXPRESSED AS
58CC                                 SQRT(AK1*AK2)/[H+])
59CC                                 AT END OF ONE ITERATION STEP
60CC      *CTO*              REAL    ACTUAL [SUM(C-12)O2], DUMMY VARIABLE
61CC      *CT1*              REAL    ACTUAL [SUM(C-12)O2], DUMMY VARIABLE
62CC      *ALKA*             REAL    ACTUAL [ALK] [EQV/L], DUMMY VARIABLE
63CC      *AKB*              REAL    1. DISSOC. CONSTANT OF BORIC ACID
64CC      *H2CO3(jpi,jpj)*     REAL    SURFACE [H2CO3] [MOLE/L]
65CC      *TWR*              REAL    2., HALF PRECISIONION
66CC      *FOURR*            REAL    4., HALF PRECISIONION
67CC      *YEAR*             REAL    NUMBER OF SECONDS IN 1 YEAR
68CC      *KRORR*         INTEGER    COUNTS ITERATIONS IN NEWTON-RAPHSON ALGO-
69CC                                 RITHM
70CC      *alpco2*            REAL    AK0, SOLUBILITY OF CO2 IN SEAWATER, DUMMY
71CC                                 VARIABLE
72CC      *BT*               REAL    TOTAL BORATE [MOLES/L]
73CC                                 ([B(OH)3]+[B(OH)4-]), DUMMY VARIABLE
74CC      *CT1*              REAL    ACTUAL VALUE OF INORG. C, DUMMY VARIABLE
75CC      *X1*               REAL    [H+] EXPRESSED AS SQRT(AK1*AK2)/[H+],
76CC                                 DUMMY VARIABLE
77CC      *PCO2*             REAL    PCO2 OF SURFACE WATER [PPM], DUMMY VARIABLE
78CC      *PA*               REAL    ATMOSPHERIC PCO2 [PPM], DUMMY VARIABLE
79CC      *FUGAC*            REAL    ENHANCEMENT FACTOR FOR GAS EXCHANGE FLUXES
80CC      *FLD*              REAL    [CO2] FLUX ATMOSPHERE -> OCEAN
81CC                                 IN [PPM/TIME STEP]
82CC      *FLU*              REAL    [CO2] FLUX OCEAN -> ATMOSPHER
83CC                                 IN [PPM/TIME STEP]
84CC      *RELW13*           REAL    RATIO [SUM((13C)O2)]/[SUM((12C)O2)]
85CC                                 IN SURFACE LAYER
86CC      *RELA13*           REAL    RATIO [SUM((13C)O2)]/[SUM((12C)O2)]
87CC                                 IN THE ATMOSPHERE
88CC      *RELW14*           REAL    RATIO [SUM((14C)O2)]/[SUM((12C)O2)]
89CC                                 IN SURFACE LAYER
90CC      *RELA14*           REAL    RATIO [SUM((14C)O2)]/[SUM((12C)O2)]
91CC                                 IN THE ATMOSPHERE
92CC
93CC   MODIFICATIONS:
94CC   --------------
95CC      original      : 1988-07 E. MAIER-REIMER      MPI HAMBURG
96CC      additions     : 1998    O. Aumont
97CC      modifications : 1999    C. Le Quere
98CC     -----------------------------------------------------------------
99CC  parameters and commons
100CC ======================
101CDIR$ NOLIST
102      USE oce_trc
103      USE trp_trc
104      USE sms
105      IMPLICIT NONE
106CDIR$ LIST
107CC----------------------------------------------------------------------
108CC local declarations
109CC ==================
110C
111      INTEGER ji, jj, it
112      REAL zexp1, zexp2
113      REAL a1, a2, a3, b2, b3, ttc, ws, krorr, alpco2
114      REAL pa, fld, flu, conveg
115      REAL oxy16, voro16, flu16
116      REAL relw13, rela13, fra13, frw13
117      REAL h,ah2,bot
118      REAL ct1,a, alka
119      REAL schmitt, caralk, bicarb
120C
121C  1. ASSIGNATION TO EXPONENTS IN THE LISS AND MERLIVAT
122C     FORMULATION OF THE GAS EXCHANGE RATE
123c -----------------------------------------------------
124C
125      zexp1 = -2./3.
126      zexp2 = -1./2. 
127      a1    = 0.17
128      a2    = 2.85
129      a3    = 5.90
130      b2    = 9.65
131      b3    = 49.3     
132C
133C  2.1 INITIALISATION OF [H+]
134C  --------------------------
135C
136      DO jj=1,jpj
137        DO ji=1,jpi
138          caralk = trn(ji,jj,1,jptal)-
139     &                 borat(ji,jj,1)/(1.+1E-8/akb3(ji,jj,1))
140          co3(ji,jj,1) = caralk-trn(ji,jj,1,jpdic)
141     &        +(1.-tmask(ji,jj,1))*.5e-3
142          bicarb = (2.*trn(ji,jj,1,jpdic)-caralk)
143          hi(ji,jj,1) = ak23(ji,jj,1)*bicarb/co3(ji,jj,1)
144        END DO
145      END DO
146     
147C
148C* 2.2 SURFACE CHEMISTRY (PCO2 AND [H+] IN
149C     SURFACE LAYER); THE RESULT OF THIS CALCULATION
150C     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
151C ---------------------------------------------------
152C
153      DO krorr = 1,15
154C
155        DO jj = 1,jpj
156          DO ji = 1,jpi
157C
158C* 2.3 TOTAL BORATE
159C -----------------
160C
161            bot = chemc(ji,jj,2)
162            ct1  = trn(ji,jj,1,jpdic)
163            alka = trn(ji,jj,1,jptal) 
164            h = hi(ji,jj,1)+(1.-tmask(ji,jj,1))*1.e-9
165            h = amax1(hi(ji,jj,1),1.E-10)
166           
167C
168C* 2.4 CALCULATE [ALK]([CO3--], [HCO3-])
169C ------------------------------------
170C
171            a=alka-
172     &          (akw3(ji,jj,1)/h-h+bot/(1.+h/akb3(ji,jj,1)))
173C
174C* 2.5 CALCULATE [H+] AND [H2CO3]
175C -----------------------------------------
176C
177            ah2=sqrt((ct1-a)**2+4*(a*ak23(ji,jj,1)/ak13(ji,jj,1))
178     &          *(2*ct1-a)) 
179            ah2=0.5*ak13(ji,jj,1)/a*((ct1-a)+ah2)
180            hi(ji,jj,1)  = ah2
181            h2co3(ji,jj) = (2*ct1-a)/(2.+ak13(ji,jj,1)/ah2)
182            conveg=((ah2-hi(ji,jj,1))/hi(ji,jj,1))**2
183     $          *tmask(ji,jj,1)
184            rconvs = rconvs+conveg
185          END DO
186        END DO
187C
188C  2.6 TEST CONVERGENCE
189C ---------------------
190C
191        IF (rconvs.LE.1.E-4) EXIT
192C
193      END DO
194C
195C
196C  2.7 CHECK ITERATION
197C  --------------------
198C
199C     IF (krorr.GE.30) THEN
200C         WRITE(numout,*) 'in h3cflx h2co3 have not converged '
201C         CALL FLUSH (numout)
202C         STOP 'h3cflx'
203C     END IF
204C
205C
206C 3. COMPUTE FLUXES
207C --------------
208     
209C
210C 3.1 FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
211C -------------------------------------------
212C
213     
214      DO jj = 1,jpj
215        DO ji = 1,jpi
216C
217          ttc = min(39.,tn(ji,jj,1))
218          schmitt= 2073.1-125.62*ttc+3.6276*ttc**2-0.043126*ttc**3
219C
220C 3.2 COMPUTE GAS EXCHANGE FOR CO2
221C --------------------------------
222C
223          if (igaswind.eq.1) then
224              ws = vatm(ji,jj)
225C
226C 3.3 this is Wanninkopf(1992) equation 8 (with chemical enhancement), 
227C     in cm/h
228C --------------------------------------------------------------------
229C
230              kgwanin(ji,jj) = (0.3*ws*ws + 2.5*(0.5246+ttc*(0.016256+
231     &            ttc*0.00049946)))*sqrt(660./schmitt)
232C
233C 3.4 CONVERT TO M/S
234C ------------------
235C
236              kgwanin(ji,jj) = kgwanin(ji,jj)/100./3600.
237C
238C 3.5 convert to mol/m2/s/uatm, alpco2(chemc(ji,jj,1)) is in 
239C      mol/L/uatm and apply ice cover
240C -----------------------------------------------------------
241C
242              kgwanin(ji,jj) = kgwanin(ji,jj)*chemc(ji,jj,1)*1.e3 * 
243     &            (1-freeze(ji,jj))
244C
245C 3.6 CLIMATOLOGICAL CASE: GAS EXCHANGE IS NOT COMPUTED BUT
246C     READ IN A FILE
247C ---------------------------------------------------------
248C
249           ELSE IF (igaswind.EQ.2) THEN
250C
251C 3.7 CONVERT TO M/S
252C ------------------
253C
254               kgwanin(ji,jj) = kgwanin(ji,jj)/raass
255     &             *(1-freeze(ji,jj))
256C
257C
258C 3.8 COUPLED RUN CASE. WIND FIELD IS COMPUTED FROM WIND STRESS
259C --------------------------------------------------------------
260C
261           ELSE IF (igaswind.EQ.3) THEN
262               
263C 3.9 LB zonal stress = 0 on east side of bassins
264C     LB uses stress from coupled runs to compute kgwanin
265C --------------------------------------------------------
266C
267               ws=sqrt((sqrt((taux(ji,jj)**2)+(tauy(ji,jj))**2))
268     $             /(1.25e-3))
269               kgwanin(ji,jj)=(0.3*(sqrt((taux(ji,jj)**2)+
270     $             (tauy(ji,jj)) **2))/(1.25e-3)  + 2.5*(0.5246 +
271     $             ttc*(0.016256 +ttc*0.00049946)))*sqrt(660./schmitt)
272C
273C 3.10 CONVERT TO M/S
274C --------------------
275C
276               kgwanin(ji,jj) = kgwanin(ji,jj)/360000.
277C
278C 3.11 convert to mol/m2/s/uatm, alpco2(chemc(ji,jj,1)) is in mol/L/uatm
279C      and apply ice cover
280C ----------------------------------------------------------------------
281C
282               kgwanin(ji,jj)=kgwanin(ji,jj)*chemc(ji,jj,1)*1.e3
283     &             *(1-freeze(ji,jj))
284           ENDIF
285         END DO
286       END DO
287C
288C
289C 3.12 COMPUTE GAS EXCHANGE COEFFICIENT FO O2 FROM LISS AND
290C      MERLIVAT EQUATIONS
291C ---------------------------------------------------------
292C
293       DO jj = 1,jpj
294         DO ji=1,jpi
295C
296           IF (igaswind.EQ.3) then
297               ws=sqrt((sqrt((taux(ji,jj)**2)+(tauy(ji,jj))**2))
298     $             /(1.25e-3))
299           ELSE
300               ws = vatm(ji,jj)
301           ENDIF
302           
303           ttc = min(39.,tn(ji,jj,1))
304           SChmitt = 1953.4-128.0*ttc+3.9918*ttc**2-0.050091*ttc**3
305C
306C 
307C
308           IF (ws.LE.3.6) THEN
309               fugaci(ji,jj) = (a1*ws)*(schmitt/660.)**zexp1
310           ENDIF
311           IF ((ws.GT.3.6).AND.(ws.LE.13.)) THEN
312               fugaci(ji,jj) = (a2*ws-b2)*(schmitt/660.)**zexp2
313           ENDIF
314           IF (ws.GT.13.) THEN
315               fugaci(ji,jj) = (a3*ws-b3)*(schmitt/660.)**zexp2
316           ENDIF
317C
318C convert to cm and apply ice cover
319C
320           fugaci(ji,jj) = fugaci(ji,jj)/100./(3600.)*
321     $         (1-freeze(ji,jj))*tmask(ji,jj,1)
322C
323#    if defined key_off_degrad
324           fugaci(ji,jj) = exp(-rfact*fugaci(ji,jj)
325     $         *facvol(ji,jj,1)/e3t(1))
326#    else
327           fugaci(ji,jj) = exp(-rfact*fugaci(ji,jj)
328     $         /e3t(1))
329#    endif
330           
331         ENDDO
332       ENDDO
333C
334       DO jj = 1,jpj
335         DO ji = 1,jpi
336C
337C Fix atmospheric CO2
338C -------------------
339C
340C for climate runs
341           pa = atcco2 
342C
343C Compute CO2 flux for the sea and air
344C ------------------------------------
345C
346           alpco2 = chemc(ji,jj,1)
347           fld = pa*tmask(ji,jj,1)*kgwanin(ji,jj)
348           flu = h2co3(ji,jj)/alpco2*tmask(ji,jj,1)*kgwanin(ji,jj)
349           tra(ji,jj,1,jpdic)= tra(ji,jj,1,jpdic)+(fld-flu)/1000./e3t(1)
350CMAFOA faire qcumul sur domaine exact
351           qcumul(1)=qcumul(1)+(fld-flu)*rfact*e1t(ji,jj)
352     $         *e2t(ji,jj)*tmask(ji,jj,1)
353C
354#    if defined key_trc_biohamocc13
355C
356C Compute C13 flux
357C ----------------
358C
359          relw13 = trn(ji,jj,1,jp13c)/
360     &       (trn(ji,jj,1,jpdic)+(1.-tmask(ji,jj,1))*1.e-20)
361          rela13 = pdb*(1.-6.4/1000.)
362          fra13  = 1.+(-10.66+0.1196*tn(ji,jj,1)-0.0003095*
363     &          tn(ji,jj,1)**2)/1000.
364          frw13  = 1.
365          tra(ji,jj,1,jp13c) = tra(ji,jj,1,jp13c)+
366     &       (fld*frw13*rela13-flu*relw13*fra13)/1000./e3t(1)
367#    endif
368C
369C Compute O2 flux 
370C ---------------
371C
372          oxy16 = trn(ji,jj,1,jpoxy)
373          voro16 = atcox
374          flu16 = (-fugaci(ji,jj)+1)*e3t(1)
375     &            *(voro16*chemc(ji,jj,3)-oxy16)*
376     &            tmask(ji,jj,1)/rfact 
377          tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy)+flu16/e3t(1)
378C
379C
380C
381C Save diagnostics
382C ----------------
383C
384#    if defined key_trc_diaadd
385          trc2d(ji,jj,1) = (fld-flu)
386          trc2d(ji,jj,2) = flu16*1000.
387          trc2d(ji,jj,3) = kgwanin(ji,jj)
388          trc2d(ji,jj,4) = (fld-flu)/(kgwanin(ji,jj)+1.E-15)
389C
390#        if defined key_trc_biohamocc13
391          trc2d(ji,jj,10) = (fld*frw13*rela13-flu*relw13*fra13)
392#        endif
393C
394#    endif
395C
396        END DO
397      END DO
398C
399#endif
400      RETURN
401      END
402
Note: See TracBrowser for help on using the repository browser.