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.
sedchem.F90 in branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/TOP_SRC/SED – NEMO

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/TOP_SRC/SED/sedchem.F90 @ 3524

Last change on this file since 3524 was 3524, checked in by gm, 11 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. add USE lib_fortran when SIGN is used (TOP,OPA,LIM2&3) ; salt flux names start with sfx_ in LIM3

  • Property svn:keywords set to Id
File size: 22.8 KB
Line 
1MODULE sedchem
2   !
3#if defined key_sed 
4   !!======================================================================
5   !!                        ***  Module sedchem  ***
6   !! sediment :   Variable for chemistry of the CO2 cycle
7   !!======================================================================
8   USE sed            ! sediment global variable
9   USE sedarr
10   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
11
12   PUBLIC sed_chem   
13
14   !! * Module variables
15   REAL(wp) :: &
16      salchl = 1./1.80655 , & ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)
17      calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
18
19   REAL(wp) :: &              ! coeff. for 1. dissoc. of carbonic acid (Millero, 1995)   
20      c10 = 3670.7   , &
21      c11 = 62.008   , &
22      c12 = 9.7944   , &
23      c13 = 0.0118   , &
24      c14 = 0.000116
25
26   REAL(wp) :: &              ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
27      c20 = 1394.7   , &
28      c21 = 4.777    , &
29      c22 = 0.       , &
30      c23 = 0.0184   , &
31      c24 = 0.000118
32
33   REAL(wp) :: &              ! coeff. for 1. dissoc. of boric acid (Dickson and Goyet, 1994)
34      cb0  = -8966.90, &
35      cb1  = -2890.53, &
36      cb2  = -77.942 , &
37      cb3  = 1.728   , &
38      cb4  = -0.0996 , &
39      cb5  = 148.0248, &
40      cb6  = 137.1942, &
41      cb7  = 1.62142 , &
42      cb8  = -24.4344, &
43      cb9  = -25.085 , &
44      cb10 = -0.2474 , &
45      cb11 = 0.053105
46
47   REAL(wp) :: &             ! borat constants
48      bor1 = 0.000232, &
49      bor2 = 1./10.811
50
51   REAL(wp) :: &             ! constants for calculate concentrations
52      st1  = 0.14     , &    ! for sulfate (Morris & Riley 1966)
53      st2  = 1./96.062, &
54      ks0  = 141.328  , &
55      ks1  = -4276.1  , &
56      ks2  = -23.093  , &
57      ks3  = -13856.  , &
58      ks4  = 324.57   , &
59      ks5  = -47.986  , &
60      ks6  = 35474.   , &
61      ks7  = -771.54  , &
62      ks8  = 114.723  , &
63      ks9  = -2698.   , &
64      ks10 = 1776.    , &
65      ks11 = 1.       , &
66      ks12 = -0.001005 
67
68   REAL(wp) :: &             ! constants for calculate concentrations
69      ft1  = 0.000067   , &  ! fluorides (Dickson & Riley 1979 )
70      ft2  = 1./18.9984 , &
71      kf0  = -12.641    , &
72      kf1  = 1590.2     , &
73      kf2  = 1.525      , &
74      kf3  = 1.0        , &
75      kf4  =-0.001005
76
77   REAL(wp) :: &             ! coeff. for dissoc. of water (Dickson and Riley, 1979 )
78      cw0 = -13847.26  , &
79      cw1 = 148.9802   , &
80      cw2 = -23.6521   , &
81      cw3 = 118.67     , &
82      cw4 = -5.977     , &
83      cw5 = 1.0495     , &
84      cw6 = -0.01615
85 
86   REAL(wp) :: &             ! coeff. for dissoc. of phosphate (Millero (1974)
87      cp10 = 115.54    , &
88      cp11 = -4576.752 , &
89      cp12 = -18.453   , &
90      cp13 = -106.736  , &
91      cp14 = 0.69171   , &
92      cp15 = -0.65643  , &
93      cp16 = -0.01844  , &
94      cp20 = 172.1033  , &
95      cp21 = -8814.715 , &
96      cp22 = -27.927   , &
97      cp23 = -160.340  , &
98      cp24 = 1.3566    , &
99      cp25 = 0.37335   , &
100      cp26 = -0.05778  , &
101      cp30 = -18.126   , &
102      cp31 = -3070.75  , &
103      cp32 = 17.27039  , &
104      cp33 = 2.81197   , &
105      cp34 = -44.99486 , &
106      cp35 = -0.09984
107
108   REAL(wp) :: &             ! coeff. for dissoc. of silicates (Millero (1974) 
109      cs10 = 117.40    , & 
110      cs11 = -8904.2   , & 
111      cs12 = -19.334   , & 
112      cs13 = -458.79   , & 
113      cs14 = 3.5913    , & 
114      cs15 = 188.74    , & 
115      cs16 = -1.5998   , & 
116      cs17 = -12.1652  , & 
117      cs18 = 0.07871   , & 
118      cs19 = 0.        , & 
119      cs20 = 1.        , & 
120      cs21 = -0.001005
121
122   REAL(wp) :: &            ! coeff. for apparent solubility equilibrium
123      ! akcc1 = -34.452 , & ! of calcite (Ingle,1975,  1800, eq. 6)
124      ! akcc2 = -39.866 , & 
125      ! akcc3 = 110.21  , & 
126      ! akcc4 = -7.5752E-6 
127      akcc1 = -171.9065 , &    ! Millero et al. 1995 from Mucci 1983
128      akcc2 = -0.077993 , & 
129      akcc3 = 2839.319  , & 
130      akcc4 = 71.595    , & 
131      akcc5 = -0.77712  , & 
132      akcc6 = 0.0028426 , & 
133      akcc7 = 178.34    , & 
134      akcc8 = -0.07711  , & 
135      akcc9 = 0.0041249
136
137   REAL(wp) :: &                 ! universal gas constant
138      rgas = 83.145              ! bar.cm3/(mol.K) or R=8.31451 J/(g.mol.K)
139
140
141   ! coeff. for seawater pressure correction (millero 95)     
142   REAL(wp), DIMENSION(8)  :: & 
143      devk1, devk2, devk3, devk4, devk5
144
145   DATA devk1/ -25.5   , -15.82    , -29.48  , -25.60    , -48.76    , -14.51   , -23.12   , -26.57   /   
146   DATA devk2/ 0.1271  , -0.0219   , 0.1622  , 0.2324    , -0.5304   , 0.1211   ,  0.1758  , 0.2020   /   
147   DATA devk3/ 0.      , 0.        , 2.608E-3, -3.6246E-3, 0.        , -0.321E-3, -2.647E-3, -3.042E-3/   
148   DATA devk4/-3.08E-3 , 1.13E-3   , -2.84E-3, -5.13E-3  , -11.76E-3 , -2.67E-3 , -5.15E-3 , -4.08E-3 /   
149   DATA devk5/0.0877E-3, -0.1475E-3, 0.      , 0.0794E-3 , -0.3692E-3, 0.0427E-3,  0.09E-3 , 0.0714E-3/
150
151
152   ! coeff. for density of sea water (Millero & Poisson 1981)
153   REAL(wp), DIMENSION(5)  :: Adsw                       
154   DATA Adsw/8.24493E-1, -4.0899E-3, 7.6438E-5 , -8.246E-7, 5.3875E-9 /
155
156   REAL(wp), DIMENSION(3)  :: Bdsw 
157   DATA Bdsw / -5.72466E-3, 1.0227E-4, -1.6546E-6 /
158
159   REAL(wp)  :: Cdsw = 4.8314E-4
160
161   REAL(wp), DIMENSION(6)  :: Ddsw                   
162   DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/
163
164CONTAINS
165
166   SUBROUTINE sed_chem( kt )
167      !!----------------------------------------------------------------------
168      !!                   ***  ROUTINE sed_chem  ***
169      !!
170      !! ** Purpose :   set chemical constants
171      !!
172      !!   History :
173      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
174      !!        !  06-04  (C. Ethe)  Re-organization
175      !!----------------------------------------------------------------------
176      !!* Arguments
177      INTEGER, INTENT(in) :: kt                     ! time step
178
179#if ! defined key_sed_off
180      INTEGER  :: ji, jj, ikt
181      REAL(wp) :: ztkel, ztc, ztc2, zpres, ztr 
182      REAL(wp) :: zsal, zsal2, zsqrt, zsal15 
183      REAL(wp) :: zis, zis2, zisqrt         
184      REAL(wp) :: zdens0, zaw, zbw, zcw   
185      REAL(wp) :: zbuf1, zbuf2 
186      REAL(wp) :: zcpexp, zcpexp2
187      REAL(wp) :: zck1p, zck2p, zck3p, zcksi   
188      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zchem_data
189
190#endif
191      !!----------------------------------------------------------------------
192
193      IF( MOD( kt - 1, nfreq ) /= 0 ) RETURN
194
195      WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt
196      WRITE(numsed,*) ' '
197
198
199#if defined key_sed_off
200      CALL sed_chem_off
201#else
202      ! reading variables
203      ALLOCATE( zchem_data(jpi,jpj,6) )   ;   zchem_data(:,:,:) = 0.
204
205      DO jj = 1,jpj
206         DO ji = 1, jpi
207            ikt = mbkt(ji,jj) 
208            IF ( tmask(ji,jj,ikt) == 1 ) THEN
209               zchem_data(ji,jj,1) = ak13 (ji,jj,ikt)
210               zchem_data(ji,jj,2) = ak23 (ji,jj,ikt)
211               zchem_data(ji,jj,3) = akb3 (ji,jj,ikt)
212               zchem_data(ji,jj,4) = akw3 (ji,jj,ikt)
213               zchem_data(ji,jj,5) = aksp (ji,jj,ikt)
214               zchem_data(ji,jj,6) = borat(ji,jj,ikt)
215            ENDIF
216         ENDDO
217      ENDDO
218
219      CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1), iarroce(1:jpoce) )
220      CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2), iarroce(1:jpoce) )
221      CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3), iarroce(1:jpoce) )
222      CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4), iarroce(1:jpoce) )
223      CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5), iarroce(1:jpoce) )
224      CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6), iarroce(1:jpoce) )
225
226      DEALLOCATE( zchem_data )
227
228      DO ji = 1, jpoce
229         ztkel   = temp(ji) + 273.16
230         ztc     = temp(ji)
231         ztc2    = ztc * ztc
232         zpres   = press(ji)
233         ! zqtt    = ztkel * 0.01
234         zsal    = salt(ji)
235         zsal2   = zsal * zsal 
236         zsqrt   = SQRT( zsal )
237         zsal15  = zsqrt * zsal
238         zlogt   = LOG( ztkel )
239         ztr     = 1./ ztkel
240         ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet)
241         zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal )
242         zis2    = zis * zis
243         zisqrt  = SQRT( zis )
244
245         ! Density of Sea Water - F(temp,sal) [kg/m3]
246         zdens0 =  Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 &
247                  + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 &
248                  + Ddsw(6) * ztc * ztc2 * ztc2
249         zaw =  Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 &
250              + Adsw(5) * ztc2 * ztc2
251         zbw =  Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2
252         zcw =  Cdsw
253         densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal
254         densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l]
255
256         ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
257         ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982)
258         ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
259         ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.)
260         ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP
261         ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286
262         ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)
263         !-----------------------------------------------------------------------------------------
264         zcpexp  = zpres / ( rgas*ztkel )
265         zcpexp2 = zpres * zcpexp
266
267         ! For Phodphates (phosphoric acid) (DOE 1994)
268         !----------------------------------------------
269         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
270            &      + ( cp15*ztr + cp16 ) * zsal
271         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
272            &      + ( cp25*ztr + cp26 ) * zsal
273         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
274            &      + ( cp34*ztr + cp35 ) * zsal
275
276         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
277         !--------------------------------------------------------
278         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
279            &      + ( cs15*ztr + cs16 ) * zis &
280            &      + ( cs17*ztr + cs18 ) * zis2 &
281            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
282
283
284         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
285         !---------------------------------------------------
286         zak1p  = EXP ( zck1p  )
287         zak2p  = EXP ( zck2p  )
288         zak3p  = EXP ( zck3p  )
289         zaksil = EXP ( zcksi  )
290
291         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
292         zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc )
293         aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
294
295         zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 )
296         zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc )
297         ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
298 
299         zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 )
300         zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc )
301         ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
302
303         zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 )
304         zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc )
305         ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
306
307         ak12s  (ji) = ak1s (ji) * ak2s (ji)
308         ak12ps (ji) = ak1ps(ji) * ak2ps(ji)
309         ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji)
310
311         calcon2(ji) = 0.01028 * ( salt(ji) / 35. ) * densSW(ji)
312      ENDDO
313
314       
315#endif
316
317   END SUBROUTINE sed_chem
318
319#if defined key_sed_off
320
321   SUBROUTINE sed_chem_off
322      !!----------------------------------------------------------------------
323      !!                   ***  ROUTINE sed_chem_off  ***
324      !!                   
325      !! ** Purpose :   compute chemical constants
326      !!
327      !!   History :
328      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
329      !!        !  06-04  (C. Ethe)  Re-organization
330      !!---------------------------------------------------------------------- 
331      !! * Local declarations
332      INTEGER :: ji
333
334      REAL(wp) ::  ztkel, ztc, ztc2, zpres, ztr 
335      REAL(wp) ::  zsal, zsal2, zsqrt, zsal15   
336      REAL(wp) ::  zis, zis2, zisqrt           
337      REAL(wp) ::  zdens0, zaw, zbw, zcw       
338      REAL(wp) ::  zchl, zst, zft, zbuf1, zbuf2 
339      REAL(wp) ::  zcpexp, zcpexp2             
340      REAL(wp) ::  zckb, zck1, zck2, zckw       
341      REAL(wp) ::  zck1p, zck2p, zck3p, zcksi   
342      REAL(wp) ::  zak1, zak2, zakb, zakw             
343      REAL(wp) ::  zaksp0, zksp, zks, zkf 
344
345
346      ! CHEMICAL CONSTANTS - DEEP OCEAN
347      !-------------------------------------
348      ! [chem constants]=mol/kg solution (or (mol/kg sol)2 for akws and aksp)
349
350      DO ji = 1, jpoce
351         ztkel   = temp(ji) + 273.16
352         ztc     = temp(ji)
353         ztc2    = ztc * ztc
354         zpres   = press(ji)
355         ! zqtt    = ztkel * 0.01
356         zsal    = salt(ji)
357         zsal2   = zsal * zsal 
358         zsqrt   = SQRT( zsal )
359         zsal15  = zsqrt * zsal
360         zlogt   = LOG( ztkel )
361         ztr     = 1./ ztkel
362         ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet)
363         zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal )
364         zis2    = zis * zis
365         zisqrt  = SQRT( zis )
366
367
368         ! Density of Sea Water - F(temp,sal) [kg/m3]
369         zdens0 =  Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 &
370                  + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 &
371                  + Ddsw(6) * ztc * ztc2 * ztc2
372         zaw =  Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 &
373              + Adsw(5) * ztc2 * ztc2
374         zbw =  Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2
375         zcw =  Cdsw
376         densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal
377         densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l]
378
379
380         ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)
381         ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982)
382         ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
383         ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.)
384         ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP
385         ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286
386         ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)
387         !-----------------------------------------------------------------------------------------
388         zcpexp  = zpres / ( rgas*ztkel )
389         zcpexp2 = zpres * zcpexp
390
391
392         ! chlorinity (WOOSTER ET AL., 1969)
393         !---------------------------------------
394         zchl = zsal * salchl
395
396         ! total sulfate concentration [mol/kg soln]
397         ! --------------------------------------
398         zst = st1 * zchl * st2
399
400         ! total fluoride concentration [mol/kg soln]
401         ! --------------------------------------
402         zft = ft1 * zchl * ft2
403
404         ! dissociation constant for carbonate (Mehrback 74 - Dickson & Millero 87)
405         !---------------------------------------------------------------------------
406         zck1 = c10*ztr - c11 + c12*zlogt - c13*zsal + c14*zsal2
407         zck2 = c20*ztr + c21 - c22*zlogt - c23*zsal + c24*zsal2
408
409         ! dissociation constant for sulfates (Dickson 1990)
410         !--------------------------------------------------
411         zks = EXP(  ks0 + ks1*ztr + ks2*zlogt &
412            &    + ( ks3*ztr + ks4 + ks5*zlogt ) * zisqrt &
413            &    + ( ks6*ztr + ks7 + ks8*zlogt ) * zis    &
414            &    +   ks9*ztr*zis*zisqrt + ks10*ztr*zis2   &
415            &    +   LOG( ks11 + ks12*zsal ) )
416
417         ! dissociation constant for fluorides (Dickson and Riley 79)
418         !--------------------------------------------------
419         zkf = EXP( kf0 + kf1*ztr + kf2*zisqrt + LOG( kf3 + kf4*zsal ) )
420
421         ! dissociation constant for borates (Doe 94)
422         !--------------------------------------------------
423         zckb = (  cb0 + cb1*zsqrt + cb2*zsal + cb3*zsal15 + cb4*zsal2) * ztr &
424            &  + ( cb5 + cb6*zsqrt + cb7*zsal) &
425            &  + ( cb8 + cb9*zsqrt + cb10*zsal) * zlogt &
426            &  +   cb11*zsqrt*ztkel + LOG( ( 1. + zst/zks + zft/zkf ) / ( 1. + zst/zks ) ) 
427
428         ! PKW (H2O) (DICKSON AND RILEY, 1979)
429         !--------------------------------------
430         zckw =   cw0*ztr + cw1 + cw2*zlogt &
431            & +( cw3*ztr + cw4 + cw5*zlogt )* zsqrt + cw6*zsal
432         
433         ! For Phodphates (phosphoric acid) (DOE 1994)
434         !----------------------------------------------
435         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
436            &      + ( cp15*ztr + cp16 ) * zsal
437         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
438            &      + ( cp25*ztr + cp26 ) * zsal
439         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
440            &      + ( cp34*ztr + cp35 ) * zsal
441
442         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
443         !--------------------------------------------------------
444         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
445            &      + ( cs15*ztr + cs16 ) * zis &
446            &      + ( cs17*ztr + cs18 ) * zis2 &
447            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
448
449         ! apparent solublity product K'SP of calcite in seawater
450         ! (S=27-43, T=2-25 DEG C) AT pres =0 (INGLE, 1975, EQ. 6)
451         ! prob: olivier a log = ln et C. Heize a LOG10(sal)
452         ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)+akcc4*tkel*tkel)
453         ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)+akcc4*tkel*tkel)
454         !--------------------------------------------------------------------
455         zaksp0 = akcc1 + akcc2*ztkel + akcc3*ztr + akcc4 * LOG10(ztkel) &
456            &  + ( akcc5 + akcc6*ztkel+ akcc7*ztr ) * zsqrt &
457            &  +  akcc8*zsal + akcc9*zsal15
458
459         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
460         !---------------------------------------------------
461         zak1   = 10**( -zck1  )
462         zak2   = 10**( -zck2  )
463         zakb   = EXP ( zckb   ) 
464         zakw   = EXP ( zckw   )
465         zksp   = 10**( zaksp0 )
466
467
468
469         ! KB of boric acid, K1,K2 of carbonic acid pressure correction
470         ! after Culberson and  AND Pytkowicz (1968) (CF. BROECKER ET AL., 1982) Millero 95
471         !--------------------------------------------------------------------------------
472         zbuf1       = - ( devk1(1) + devk2(1)*ztc + devk3(1)*ztc2 )
473         zbuf2       = 0.5 * ( devk4(1) + devk5(1)*ztc )
474         ak1s(ji)    = zak1 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
475
476         zbuf1       = -( devk1(2) + devk2(2)*ztc + devk3(2)*ztc2 )
477         zbuf2       = 0.5 * ( devk4(2) + devk5(2)*ztc )
478         ak2s(ji)    = zak2 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
479
480         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
481         zbuf2       = 0.5 * ( devk4(3) + devk5(3) * ztc )
482         akbs(ji)    = zakb * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
483
484         zbuf1       = - ( devk1(4) + devk2(4)*ztc + devk3(4)*ztc2 )
485         zbuf2       = 0.5 * ( devk4(4) + devk5(4)*ztc )
486         akws(ji)    = zakw * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
487
488
489         ! APPARENT SOLUBILITY PRODUCT K''SP OF CALCITE (OR ARAGONITE)
490         ! AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970)
491         ! (P. 1285) AND BERNER (1976)
492         !-----------------------------------------------------------------
493         ! aksp(ji) = aksp0*exp(zcpexp*(devks-devkst*tc))
494         ! or following Mucci
495         zbuf1      = - ( devk1(5) + devk2(5)*ztc + devk3(5)*ztc2 )
496         zbuf2      = 0.5 *( devk4(5) + devk5(5)*ztc )
497         aksps(ji)   = zksp * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
498
499         ! For Phodphates (phosphoric acid) (DOE 1994)
500         !----------------------------------------------
501         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
502            &      + ( cp15*ztr + cp16 ) * zsal
503         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
504            &      + ( cp25*ztr + cp26 ) * zsal
505         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
506            &      + ( cp34*ztr + cp35 ) * zsal
507
508         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
509         !--------------------------------------------------------
510         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
511            &      + ( cs15*ztr + cs16 ) * zis &
512            &      + ( cs17*ztr + cs18 ) * zis2 &
513            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
514
515
516         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
517         !---------------------------------------------------
518         zak1p  = EXP ( zck1p  )
519         zak2p  = EXP ( zck2p  )
520         zak3p  = EXP ( zck3p  )
521         zaksil = EXP ( zcksi  )
522
523         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
524         zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc )
525         aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
526
527         zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 )
528         zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc )
529         ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
530 
531         zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 )
532         zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc )
533         ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
534
535         zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 )
536         zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc )
537         ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
538
539         ! total borat concentration. [mol/l]
540         ! or from Millero 1995 [mol/l] : borat(l) = 0.000416_8*(sal/35._8)*densSW(l)
541         ! --------------------------------------------------------------------------
542         borats(ji) = bor1 * zchl * bor2 * densSW(ji)
543
544         ak12s  (ji) = ak1s (ji) * ak2s (ji)
545         ak12ps (ji) = ak1ps(ji) * ak2ps(ji)
546         ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji)
547
548         calcon2(ji) = 0.01028 * ( zsal / 35. ) * densSW(ji)
549
550      ENDDO
551
552   END SUBROUTINE sed_chem_off
553
554#endif
555
556#else
557   !!======================================================================
558   !! MODULE sedchem  :   Dummy module
559   !!======================================================================
560CONTAINS
561   SUBROUTINE sed_chem( kt )         ! Empty routine
562      INTEGER, INTENT(in) :: kt
563      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
564   END SUBROUTINE sed_chem
565
566   !!======================================================================
567
568#endif
569
570END MODULE sedchem
Note: See TracBrowser for help on using the repository browser.