source: branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90 @ 775

Last change on this file since 775 was 775, checked in by gm, 13 years ago

dev_001_GM - PISCES in F90 : encapsulation of all p4z…F files in module F90 + doctor norme for local variables - compilation OK

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.7 KB
Line 
1MODULE p4zche
2   !!======================================================================
3   !!                         ***  MODULE p4zche  ***
4   !! TOP :   PISCES Sea water chemistry computed following OCMIP protocol
5   !!======================================================================
6   !! History :    -   !  1988     (E. Maier-Reimer)  Original code
7   !!              -   !  1998     (O. Aumont)  addition
8   !!              -   !  1999     (C. Le Quere)  modification
9   !!             1.0  !  2004     (O. Aumont)  modification
10   !!              -   !  2006     (R. Gangsto)  modification
11   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
12   !!----------------------------------------------------------------------
13#if defined key_pisces
14   !!----------------------------------------------------------------------
15   !!   'key_pisces'                                       PISCES bio-model
16   !!----------------------------------------------------------------------
17   !!   p4z_che        :  Sea water chemistry computed following OCMIP protocol
18   !!----------------------------------------------------------------------
19   USE oce_trc         !
20   USE trp_trc         !
21   USE sms             !
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   p4z_che    ! called in p4zprg.F90
27
28   !!* Substitution
29#include "domzgr_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
32   !! $Header:$
33   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38   SUBROUTINE p4z_che
39      !!---------------------------------------------------------------------
40      !!                     ***  ROUTINE p4z_che  ***
41      !!
42      !! ** Purpose :   Sea water chemistry computed following OCMIP protocol
43      !!
44      !! ** Method  : - ...
45      !!---------------------------------------------------------------------
46      INTEGER  ::   ji, jj, jk
47      REAL(wp) ::   ztkel, zsal , zqtt  , zbuf1 , zbuf2
48      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zcek0, zoxy  , zcpexp2
49      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
50      REAL(wp) ::   zlqtt, zqtt2, zsal15, zis   , zis2 , zisqrt
51      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
52      REAL(wp) ::   zckp1, zckp2, zckp3 , zcksi , zakp1, zakp2 , zakp3, zaksi
53      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaks , zakf  , zaksp1
54      !!---------------------------------------------------------------------
55
56      ! CHEMICAL CONSTANTS - SURFACE LAYER
57      ! ----------------------------------
58
59      DO jj = 1, jpj
60         DO ji = 1, jpi
61
62            !                             ! SET ABSOLUTE TEMPERATURE
63            ztkel = tn(ji,jj,1) + 273.16
64            zqtt  = ztkel * 0.01
65            zqtt2 = zqtt * zqtt
66            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
67            zlqtt = LOG( zqtt )
68
69            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
70            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
71            zcek0 = c00 + c01 / zqtt + c02 * zlqtt + zsal * ( c03 + c04 * zqtt + c05 * zqtt2 )
72            zcek1 = ca0 + ca1 / zqtt + ca2 * zlqtt + ca3 * zqtt2 + zsal*( ca4 + ca5 * zqtt + ca6 * zqtt2 )
73
74            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970)
75            zoxy  = ox0 + ox1 / zqtt + ox2 * zlqtt + zsal * ( ox3 + ox4 * zqtt + ox5 * zqtt2 )
76
77            !                             ! SET SOLUBILITIES OF O2 AND CO2
78            chemc(ji,jj,1) = EXP( zcek0 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
79            chemc(ji,jj,2) = EXP( zoxy  ) * oxyco
80            chemc(ji,jj,3) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.
81
82         END DO
83      END DO
84
85      ! CHEMICAL CONSTANTS - DEEP OCEAN
86      ! -------------------------------
87
88      DO jk = 1, jpk
89         DO jj = 1, jpj
90            DO ji = 1, jpi
91
92               ! SET PRESSION
93               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
94
95               ! SET ABSOLUTE TEMPERATURE
96               ztkel   = tn(ji,jj,jk) + 273.16
97               zqtt    = ztkel * 0.01
98               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
99               zsqrt  = SQRT( zsal )
100               zsal15  = zsqrt * zsal
101               zlogt  = LOG( ztkel )
102               ztr    = 1. / ztkel
103               zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
104               zis2   = zis * zis
105               zisqrt = SQRT( zis )
106               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
107
108               ! CHLORINITY (WOOSTER ET AL., 1969)
109               zcl     = zsal * salchl
110
111               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
112               zst     = st1 * zcl * st2
113
114               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
115               zft     = ft1 * zcl * ft2
116
117               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
118               zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           &
119                  &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   &
120                  &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      &
121                  &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  )
122
123               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
124               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  )
125
126               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
127               zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   &
128                  &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            &
129                  &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             &
130                  &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  )
131!!gm zsal**2 to be replaced by a *...
132               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal**2
133               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
134
135               ! PKW (H2O) (DICKSON AND RILEY, 1979)
136               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
137
138               ! DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale)
139               zckp1 = cp10 + cp11 * ztr + cp12 * zlogt + zsqrt * ( cp13 * ztr + cp14 ) + zsal * ( cp15 * ztr + cp16 )
140               zckp2 = cp20 + cp21 * ztr + cp22 * zlogt + zsqrt * ( cp23 * ztr + cp24 ) + zsal * ( cp25 * ztr + cp26 )
141               zckp3 = cp30 + cp31 * ztr                + zsqrt * ( cp32 * ztr + cp33 ) + zsal * ( cp34 * ztr + cp35 )
142               zcksi = cs10 + cs11 * ztr + cs12 * zlogt + zisqrt* ( cs13 * ztr + cs14 ) + zis * ( cs15 * ztr + cs16 )  &
143                  &                                    + zis2  * ( cs17 * ztr + cs18 ) + LOG( 1.   + cs19 * zsal )     &
144                  &                                                                    + LOG( cs20 + cs21 * zsal )
145
146               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
147               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
148               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
149                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
150
151               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
152               zak1    = 10**(zck1)
153               zak2    = 10**(zck2)
154               zakb    = EXP( zckb  )
155               zakp1   = EXP( zckp1 )
156               zakp2   = EXP( zckp2 )
157               zakp3   = EXP( zckp3 )
158               zaksi   = EXP( zcksi )
159               zakw    = EXP( zckw )
160               zaksp1  = 10**(zaksp0)
161               zaks    = exp( zcks )
162               zakf    = exp( zckf )
163
164               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
165               !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
166               !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
167               !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
168               !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
169               !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
170               !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
171               !        & GIESKES (1970), P. 1285-1286 (THE SMALL
172               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
173               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
174               zcpexp  = zpres /(rgas*ztkel)
175               zcpexp2 = zpres * zpres/(rgas*ztkel)
176
177               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
178               !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
179               !        (CF. BROECKER ET AL., 1982)
180               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
181               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
182               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
183
184               zbuf1  = -(devk1(1)+devk2(1)*ztc+devk3(1)*ztc*ztc)
185               zbuf2  = 0.5*(devk4(1)+devk5(1)*ztc)
186               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
187
188               zbuf1  =     - ( devk1(2) + devk2(2) * ztc + devk3(2) * ztc * ztc )
189               zbuf2  = 0.5 * ( devk4(2) + devk5(2) * ztc )
190               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
191
192               zbuf1  =     - ( devk1(4) + devk2(4) * ztc + devk3(4) * ztc * ztc )
193               zbuf2  = 0.5 * ( devk4(4) + devk5(4) * ztc )
194               akp13(ji,jj,jk) = zakp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
195
196               zbuf1  =     - ( devk1(5) + devk2(5) * ztc + devk3(5) * ztc * ztc )
197               zbuf2  = 0.5 * ( devk4(5) + devk5(5) * ztc )
198               akp23(ji,jj,jk) = zakp2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
199
200               zbuf1  =     - ( devk1(6) + devk2(6) * ztc + devk3(6) * ztc * ztc )
201               zbuf2  = 0.5 * ( devk4(6) + devk5(6) * ztc )
202               akp33(ji,jj,jk) = zakp3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
203
204               zbuf1  =     - ( devk1(7) + devk2(7) * ztc + devk3(7) * ztc * ztc )
205               zbuf2  = 0.5 * ( devk4(7) + devk5(7) * ztc )
206               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
207
208               !  Ksi
209               !            aksi3(ji,jj,jk) = zaksi
210               !
211               !  Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys)
212               !  "deltaVsi and deltaKsi have been estimated from the value of boric acid"
213               zbuf1  =     - ( devk1(3) + devk2(3) * ztc + devk3(3) * ztc * ztc )
214               zbuf2  = 0.5 * ( devk4(3) + devk5(3) * ztc )
215               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
216
217               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE
218               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
219               !        (P. 1285) AND BERNER (1976)
220               zbuf1  =     - ( devk1(8) + devk2(8) * ztc + devk3(8) * ztc * ztc )
221               zbuf2  = 0.5 * ( devk4(8) + devk5(8) * ztc )
222               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
223
224               !  Pressure correction for sulfate and fluoride
225               zbuf1  =     - ( devk1(9) + devk2(9) * ztc + devk3(9) * ztc * ztc )
226               zbuf2  = 0.5 * ( devk4(9) + devk5(9) * ztc )
227               aks3(ji,jj,jk) = zaks   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
228
229               zbuf1  =     - ( devk1(10) + devk2(10) * ztc + devk3(10) * ztc * ztc )
230               zbuf2  = 0.5 * ( devk4(10) + devk5(10) * ztc )
231               akf3(ji,jj,jk) = zakf   * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
232
233               ! TOTAL BORATE CONCENTR. [MOLES/L]
234               borat(ji,jj,jk) = bor1 * zcl * bor2
235
236               ! Iron and SIO3 saturation concentration from ...
237               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
238               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) )
239
240            END DO
241         END DO
242      END DO
243      !
244   END SUBROUTINE p4z_che
245
246#else
247   !!======================================================================
248   !!  Dummy module :                                   No PISCES bio-model
249   !!======================================================================
250CONTAINS
251   SUBROUTINE p4z_che( kt )                   ! Empty routine
252      INTEGER, INTENT( in ) ::   kt
253      WRITE(*,*) 'p4z_che: You should not have seen this print! error?', kt
254   END SUBROUTINE p4z_che
255#endif 
256
257   !!======================================================================
258END MODULE  p4zche
Note: See TracBrowser for help on using the repository browser.