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 @ 186

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

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