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 tags/nemo_v1_04/NEMO/TOP_SRC/SMS – NEMO

source: tags/nemo_v1_04/NEMO/TOP_SRC/SMS/h3cflx.F @ 280

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

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