Changeset 775 for branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90
- Timestamp:
- 2007-12-19T14:45:15+01:00 (16 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/TOP_SRC/PISCES_SMS/p4zche.F90
r774 r775 1 2 CCC $Header$ 3 CCC TOP 1.0 , LOCEAN-IPSL (2005) 4 C This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 5 C --------------------------------------------------------------------------- 6 CDIR$ LIST 7 SUBROUTINE p4zche 8 #if defined key_top && defined key_pisces 9 CCC--------------------------------------------------------------------- 10 CCC 11 CCC ROUTINE p4zche : PISCES MODEL 12 CCC ***************************** 13 CCC 14 CCC PURPOSE. 15 CCC -------- 16 CCC *P4ZCHE* : Sea water chemistry computed following OCMIP protocol 17 CCC 18 CCC 19 CC EXTERNALS. 20 CC ---------- 21 CC rhop 22 CC 23 CC MODIFICATIONS: 24 CC -------------- 25 CC original : 1988 E. Maier-Reimer 26 CC additions : 1998 O. Aumont 27 CC modifications : 1999 C. Le Quere 28 CC modifications : 2004 O. Aumont 29 CC modifications : 2006 R. Gangsto 30 CC---------------------------------------------------------------------- 31 CC parameters and commons 32 CC ====================== 33 CDIR$ nolist 34 USE oce_trc 35 USE trp_trc 36 USE sms 37 IMPLICIT NONE 1 MODULE 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 38 29 #include "domzgr_substitute.h90" 39 CDIR$ list 40 CC---------------------------------------------------------------------- 41 CC local declarations 42 CC ================== 43 C 44 INTEGER ji, jj, jk 45 REAL tkel, sal, qtt, zbuf1, zbuf2 46 REAL pres, tc, cl, cpexp, cek0, oxy, cpexp2 47 REAL zsqrt, ztr, zlogt, cek1 48 REAL zqtt, qtt2, sal15, zis, zis2, zisqrt 49 REAL ckb, ck1, ck2, ckw, ak1, ak2, akb, aksp0, akw 50 REAL ckp1, ckp2, ckp3, cksi, akp1, akp2, akp3, aksi 51 REAL st, ft, cks, ckf, aks, akf, aksp1 52 53 C 54 C* 1. CHEMICAL CONSTANTS - SURFACE LAYER 55 C --------------------------------------- 56 C 57 DO jj = 1,jpj 58 DO ji = 1,jpi 59 C 60 C* 1.1 SET ABSOLUTE TEMPERATURE 61 C ------------------------------ 62 C 63 tkel = tn(ji,jj,1)+273.16 64 qtt = tkel*0.01 65 qtt2=qtt*qtt 66 sal = sn(ji,jj,1) + (1.-tmask(ji,jj,1))*35. 67 zqtt=log(qtt) 68 C 69 C* 1.2 LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 70 C AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 71 C ------------------------------------------------------- 72 C 73 cek0 = c00+c01/qtt+c02*zqtt+sal*(c03+c04*qtt+c05*qtt2) 74 cek1 = ca0+ca1/qtt+ca2*zqtt+ca3*qtt2+sal*(ca4 75 & +ca5*qtt+ca6*qtt2) 76 C 77 C* 1.3 LN(K0) OF SOLUBILITY OF O2 and N2 (EQ. 4, WEISS, 1970) 78 C ------------------------------------------------------------ 79 C 80 oxy = ox0+ox1/qtt+ox2*zqtt+sal*(ox3+ox4*qtt+ox5*qtt2) 81 C 82 C* 1.4 SET SOLUBILITIES OF O2 AND CO2 83 C ----------------------------------- 84 C 85 chemc(ji,jj,1) = exp(cek0)*1.E-6*rhop(ji,jj,1)/1000. 86 chemc(ji,jj,2) = exp(oxy)*oxyco 87 chemc(ji,jj,3) = exp(cek1)*1.E-6*rhop(ji,jj,1)/1000. 88 C 89 ENDDO 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 36 CONTAINS 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 90 83 END DO 91 C 92 C* 2 CHEMICAL CONSTANTS - DEEP OCEAN 93 C ------------------------------------- 94 C 95 DO jk = 1,jpk 96 DO jj = 1,jpj 97 DO ji = 1,jpi 98 C 99 C* 2.1 SET PRESSION 100 C ----------------- 101 C 102 pres = 1.025e-1*fsdept(ji,jj,jk) 103 C 104 C* 2.2 SET ABSOLUTE TEMPERATURE 105 C ------------------------------ 106 C 107 tkel = tn(ji,jj,jk)+273.16 108 qtt = tkel*0.01 109 sal = sn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*35. 110 zsqrt = sqrt(sal) 111 sal15 = zsqrt*sal 112 zlogt = log(tkel) 113 ztr = 1./tkel 114 zis = 19.924*sal/(1000.-1.005*sal) 115 zis2 = zis*zis 116 zisqrt = sqrt(zis) 117 tc = tn(ji,jj,jk) + (1.-tmask(ji,jj,jk))*20. 118 C 119 C* 2.3 CHLORINITY (WOOSTER ET AL., 1969) 120 C --------------------------------------- 121 C 122 cl = sal*salchl 123 C 124 C* 2.4 TOTAL SULFATE CONCENTR. [MOLES/kg soln] 125 C -------------------------------------------- 126 C 127 st = st1*cl*st2 128 C 129 C* 2.5 TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 130 C --------------------------------------------- 131 C 132 ft = ft1*cl*ft2 133 C 134 C* 2.6 DISSOCIATION CONSTANT FOR SULFATES 135 C on free H scale (Dickson 1990) 136 C ------------------------------------------------------- 137 C 138 cks=exp(ks1*ztr+ks0+ks2*zlogt+(ks3*ztr+ks4+ks5*zlogt) 139 & *zisqrt+(ks6*ztr+ks7+ks8*zlogt)*zis+ks9*ztr*zis 140 & *zisqrt+ks10*ztr*zis2+log(ks11+ks12*sal)) 141 C 142 C* 2.7 DISSOCIATION CONSTANT FOR FLUORIDES 143 C on free H scale (Dickson and Riley 79) 144 C ------------------------------------------------------- 145 C 146 ckf=exp(kf1*ztr+kf0+kf2*zisqrt+log(kf3+kf4*sal)) 147 148 C 149 C* 2.4 DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 150 C ------------------------------------------------------- 151 C 152 ckb = (cb0+cb1*zsqrt+cb2*sal+cb3*sal15+cb4*sal*sal)*ztr 153 & +(cb5+cb6*zsqrt+cb7*sal)+ 154 & (cb8+cb9*zsqrt+cb10*sal)*zlogt+cb11*zsqrt*tkel 155 & +log((1.+st/cks+ft/ckf)/(1.+st/cks)) 156 ck1 = c10*ztr+c11+c12*zlogt+c13*sal+c14*sal**2 157 ck2 = c20*ztr+c21+c22*sal+c23*sal**2 158 C 159 C* 2.5 PKW (H2O) (DICKSON AND RILEY, 1979) 160 C ----------------------------------------- 161 C 162 ckw = cw0*ztr+cw1+cw2*zlogt+(cw3*ztr+cw4+cw5*zlogt)* 163 & zsqrt+cw6*sal 164 165 C 166 C 167 C* 2.10 DISSOCIATION CONSTANT FOR PHOSPHATE AND SILICATE (seawater scale) 168 C --------------------------------------------------------------------- 169 C 170 ckp1 = cp10+cp11*ztr+cp12*zlogt+zsqrt*(cp13*ztr 171 & +cp14)+sal*(cp15*ztr+cp16) 172 ckp2 = cp20+cp21*ztr+cp22*zlogt+zsqrt*(cp23*ztr 173 & +cp24)+sal*(cp25*ztr+cp26) 174 ckp3 = cp30+cp31*ztr+zsqrt*(cp32*ztr 175 & +cp33)+sal*(cp34*ztr+cp35) 176 cksi = cs10+cs11*ztr+cs12*zlogt+zisqrt*(cs13*ztr 177 & +cs14)+zis*(cs15*ztr+cs16)+zis2*(cs17*ztr 178 & +cs18)+log(1.+cs19*sal) 179 & +log(cs20+cs21*sal) 180 181 C 182 C*2.7 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 183 C (S=27-43, T=2-25 DEG C) AT pres =0 (ATMOSPH. PRESSURE) 184 C (MUCCI 1983) 185 C ------------------------------------------------------------- 186 C 187 aksp0 = akcc1+akcc2*tkel+akcc3*ztr+akcc4*log10(tkel)+ 188 & (akcc5+akcc6*tkel+ 189 & akcc7*ztr)*zsqrt+akcc8*sal+akcc9*sal15 190 191 192 C 193 C* 2.6 K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 194 C ----------------------------------------------------------------- 195 C 196 ak1 = 10**(ck1) 197 ak2 = 10**(ck2) 198 akb = exp(ckb) 199 akp1 = exp(ckp1) 200 akp2 = exp(ckp2) 201 akp3 = exp(ckp3) 202 aksi = exp(cksi) 203 akw = exp(ckw) 204 aksp1 = 10**(aksp0) 205 aks = exp(cks) 206 akf = exp(ckf) 207 208 209 C 210 C* 2.8 FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970) 211 C (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE 212 C IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 213 C TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN 214 C DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 215 C MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 216 C WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 217 C AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL 218 C FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 219 C SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON 220 C P. 1285)) 221 C ----------------------------------------------------------- 222 C 223 cpexp = pres /(rgas*tkel) 224 cpexp2 = pres * pres/(rgas*tkel) 225 C 226 C* 2.9 KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 227 C CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 228 C (CF. BROECKER ET AL., 1982) 229 C -------------------------------------------------------- 230 C 231 zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc) 232 zbuf2 = 0.5*(devk4(3)+devk5(3)*tc) 233 akb3(ji,jj,jk) = akb*exp(zbuf1*cpexp+zbuf2*cpexp2) 234 235 zbuf1 = -(devk1(1)+devk2(1)*tc+devk3(1)*tc*tc) 236 zbuf2 = 0.5*(devk4(1)+devk5(1)*tc) 237 ak13(ji,jj,jk) = ak1*exp(zbuf1*cpexp+zbuf2*cpexp2) 238 239 zbuf1 = -(devk1(2)+devk2(2)*tc+devk3(2)*tc*tc) 240 zbuf2 = 0.5*(devk4(2)+devk5(2)*tc) 241 ak23(ji,jj,jk) = ak2*exp(zbuf1*cpexp+zbuf2*cpexp2) 242 243 zbuf1 = -(devk1(4)+devk2(4)*tc+devk3(4)*tc*tc) 244 zbuf2 = 0.5*(devk4(4)+devk5(4)*tc) 245 akp13(ji,jj,jk) = akp1*exp(zbuf1*cpexp+zbuf2*cpexp2) 246 247 zbuf1 = -(devk1(5)+devk2(5)*tc+devk3(5)*tc*tc) 248 zbuf2 = 0.5*(devk4(5)+devk5(5)*tc) 249 akp23(ji,jj,jk) = akp2*exp(zbuf1*cpexp+zbuf2*cpexp2) 250 251 zbuf1 = -(devk1(6)+devk2(6)*tc+devk3(6)*tc*tc) 252 zbuf2 = 0.5*(devk4(6)+devk5(6)*tc) 253 akp33(ji,jj,jk) = akp3*exp(zbuf1*cpexp+zbuf2*cpexp2) 254 255 zbuf1 = -(devk1(7)+devk2(7)*tc+devk3(7)*tc*tc) 256 zbuf2 = 0.5*(devk4(7)+devk5(7)*tc) 257 akw3(ji,jj,jk) = akw*exp(zbuf1*cpexp+zbuf2*cpexp2) 258 259 C Ksi 260 C aksi3(ji,jj,jk) = aksi 261 C 262 C Or using coefficient of borates (cf millero 95+ corrected version html doc co2sys) 263 C "deltaVsi and deltaKsi have been estimated from the value of boric acid" 264 C 265 zbuf1 = -(devk1(3)+devk2(3)*tc+devk3(3)*tc*tc) 266 zbuf2 = 0.5*(devk4(3)+devk5(3)*tc) 267 aksi3(ji,jj,jk) = aksi*exp(zbuf1*cpexp+zbuf2*cpexp2) 268 269 C 270 C 271 C* 2.15 APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 272 C AS FUNCTION OF PRESSURE FOLLOWING MILLERO 273 C (P. 1285) AND BERNER (1976) 274 C ------------------------------------------------- 275 276 zbuf1 = -(devk1(8)+devk2(8)*tc+devk3(8)*tc*tc) 277 zbuf2 = 0.5*(devk4(8)+devk5(8)*tc) 278 aksp(ji,jj,jk) = aksp1*exp(zbuf1*cpexp+zbuf2*cpexp2) 279 280 C Pressure correction for sulfate and fluoride 281 C 282 zbuf1 = -(devk1(9)+devk2(9)*tc+devk3(9)*tc*tc) 283 zbuf2 = 0.5*(devk4(9)+devk5(9)*tc) 284 aks3(ji,jj,jk) = aks*exp(zbuf1*cpexp+zbuf2*cpexp2) 285 286 zbuf1 = -(devk1(10)+devk2(10)*tc+devk3(10)*tc*tc) 287 zbuf2 = 0.5*(devk4(10)+devk5(10)*tc) 288 akf3(ji,jj,jk) = akf*exp(zbuf1*cpexp+zbuf2*cpexp2) 289 290 291 C 292 C* 2.11 TOTAL BORATE CONCENTR. [MOLES/L] 293 C -------------------------------------- 294 C 295 borat(ji,jj,jk) = bor1*cl*bor2 296 C 297 C 2.12 Iron and SIO3 saturation concentration from ... 298 C ---------------------------------------------------- 299 C 300 sio3eq(ji,jj,jk)=exp(log(10.)*(6.44-968./tkel))*1E-6 301 fekeq(ji,jj,jk)=10**(17.27-1565.7/(273.15+tc)) 302 C 303 ENDDO 304 ENDDO 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 305 242 END DO 306 C 307 #endif 308 C 309 RETURN 310 END 243 ! 244 END SUBROUTINE p4z_che 245 246 #else 247 !!====================================================================== 248 !! Dummy module : No PISCES bio-model 249 !!====================================================================== 250 CONTAINS 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 !!====================================================================== 258 END MODULE p4zche
Note: See TracChangeset
for help on using the changeset viewer.