- Timestamp:
- 2018-11-16T16:16:27+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/dev_r9950_old_tidal_mixing/src/TOP/PISCES/SED/sedchem.F90
r9950 r10324 1 1 MODULE sedchem 2 2 3 #if defined key_sed4 3 !!====================================================================== 5 4 !! *** Module sedchem *** … … 9 8 USE sed ! sediment global variable 10 9 USE sedarr 10 USE eosbn2, ONLY : neos 11 USE lib_mpp ! distribued memory computing library 12 13 IMPLICIT NONE 14 PRIVATE 11 15 12 16 !! * Accessibility 13 PUBLIC sed_chem 17 PUBLIC sed_chem 18 PUBLIC ahini_for_at_sed ! 19 PUBLIC solve_at_general_sed ! 20 21 ! Maximum number of iterations for each method 22 INTEGER, PARAMETER :: jp_maxniter_atgen = 20 23 REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 14 24 15 25 !! * Module variables 16 26 REAL(wp) :: & 17 salchl = 1./1.80655 , & ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)18 27 calcon = 1.03E-2 ! mean calcite concentration [Ca2+] in sea water [mole/kg solution] 19 28 20 REAL(wp) :: & ! coeff. for 1. dissoc. of carbonic acid (Millero, 1995) 21 c10 = 3670.7 , & 22 c11 = 62.008 , & 23 c12 = 9.7944 , & 24 c13 = 0.0118 , & 25 c14 = 0.000116 26 27 REAL(wp) :: & ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995) 28 c20 = 1394.7 , & 29 c21 = 4.777 , & 30 c22 = 0. , & 31 c23 = 0.0184 , & 32 c24 = 0.000118 33 34 REAL(wp) :: & ! coeff. for 1. dissoc. of boric acid (Dickson and Goyet, 1994) 35 cb0 = -8966.90, & 36 cb1 = -2890.53, & 37 cb2 = -77.942 , & 38 cb3 = 1.728 , & 39 cb4 = -0.0996 , & 40 cb5 = 148.0248, & 41 cb6 = 137.1942, & 42 cb7 = 1.62142 , & 43 cb8 = -24.4344, & 44 cb9 = -25.085 , & 45 cb10 = -0.2474 , & 46 cb11 = 0.053105 47 48 REAL(wp) :: & ! borat constants 49 bor1 = 0.000232, & 50 bor2 = 1./10.811 51 52 REAL(wp) :: & ! constants for calculate concentrations 53 st1 = 0.14 , & ! for sulfate (Morris & Riley 1966) 54 st2 = 1./96.062, & 55 ks0 = 141.328 , & 56 ks1 = -4276.1 , & 57 ks2 = -23.093 , & 58 ks3 = -13856. , & 59 ks4 = 324.57 , & 60 ks5 = -47.986 , & 61 ks6 = 35474. , & 62 ks7 = -771.54 , & 63 ks8 = 114.723 , & 64 ks9 = -2698. , & 65 ks10 = 1776. , & 66 ks11 = 1. , & 67 ks12 = -0.001005 68 69 REAL(wp) :: & ! constants for calculate concentrations 70 ft1 = 0.000067 , & ! fluorides (Dickson & Riley 1979 ) 71 ft2 = 1./18.9984 , & 72 kf0 = -12.641 , & 73 kf1 = 1590.2 , & 74 kf2 = 1.525 , & 75 kf3 = 1.0 , & 76 kf4 =-0.001005 77 78 REAL(wp) :: & ! coeff. for dissoc. of water (Dickson and Riley, 1979 ) 79 cw0 = -13847.26 , & 80 cw1 = 148.9802 , & 81 cw2 = -23.6521 , & 82 cw3 = 118.67 , & 83 cw4 = -5.977 , & 84 cw5 = 1.0495 , & 85 cw6 = -0.01615 86 87 REAL(wp) :: & ! coeff. for dissoc. of phosphate (Millero (1974) 88 cp10 = 115.54 , & 89 cp11 = -4576.752 , & 90 cp12 = -18.453 , & 91 cp13 = -106.736 , & 92 cp14 = 0.69171 , & 93 cp15 = -0.65643 , & 94 cp16 = -0.01844 , & 95 cp20 = 172.1033 , & 96 cp21 = -8814.715 , & 97 cp22 = -27.927 , & 98 cp23 = -160.340 , & 99 cp24 = 1.3566 , & 100 cp25 = 0.37335 , & 101 cp26 = -0.05778 , & 102 cp30 = -18.126 , & 103 cp31 = -3070.75 , & 104 cp32 = 17.27039 , & 105 cp33 = 2.81197 , & 106 cp34 = -44.99486 , & 107 cp35 = -0.09984 108 109 REAL(wp) :: & ! coeff. for dissoc. of silicates (Millero (1974) 110 cs10 = 117.40 , & 111 cs11 = -8904.2 , & 112 cs12 = -19.334 , & 113 cs13 = -458.79 , & 114 cs14 = 3.5913 , & 115 cs15 = 188.74 , & 116 cs16 = -1.5998 , & 117 cs17 = -12.1652 , & 118 cs18 = 0.07871 , & 119 cs19 = 0. , & 120 cs20 = 1. , & 121 cs21 = -0.001005 122 123 REAL(wp) :: & ! coeff. for apparent solubility equilibrium 124 ! akcc1 = -34.452 , & ! of calcite (Ingle,1975, 1800, eq. 6) 125 ! akcc2 = -39.866 , & 126 ! akcc3 = 110.21 , & 127 ! akcc4 = -7.5752E-6 128 akcc1 = -171.9065 , & ! Millero et al. 1995 from Mucci 1983 129 akcc2 = -0.077993 , & 130 akcc3 = 2839.319 , & 131 akcc4 = 71.595 , & 132 akcc5 = -0.77712 , & 133 akcc6 = 0.0028426 , & 134 akcc7 = 178.34 , & 135 akcc8 = -0.07711 , & 136 akcc9 = 0.0041249 137 138 REAL(wp) :: & ! universal gas constant 139 rgas = 83.145 ! bar.cm3/(mol.K) or R=8.31451 J/(g.mol.K) 140 141 142 ! coeff. for seawater pressure correction (millero 95) 143 REAL(wp), DIMENSION(8) :: & 144 devk1, devk2, devk3, devk4, devk5 145 146 DATA devk1/ -25.5 , -15.82 , -29.48 , -25.60 , -48.76 , -14.51 , -23.12 , -26.57 / 147 DATA devk2/ 0.1271 , -0.0219 , 0.1622 , 0.2324 , -0.5304 , 0.1211 , 0.1758 , 0.2020 / 148 DATA devk3/ 0. , 0. , 2.608E-3, -3.6246E-3, 0. , -0.321E-3, -2.647E-3, -3.042E-3/ 149 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 / 150 DATA devk5/0.0877E-3, -0.1475E-3, 0. , 0.0794E-3 , -0.3692E-3, 0.0427E-3, 0.09E-3 , 0.0714E-3/ 151 29 REAL(wp) :: rgas = 83.14472 ! universal gas constants 152 30 153 31 ! coeff. for density of sea water (Millero & Poisson 1981) … … 162 40 REAL(wp), DIMENSION(6) :: Ddsw 163 41 DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/ 42 43 REAL(wp) :: devk10 = -25.5 44 REAL(wp) :: devk11 = -15.82 45 REAL(wp) :: devk12 = -29.48 46 REAL(wp) :: devk13 = -20.02 47 REAL(wp) :: devk14 = -18.03 48 REAL(wp) :: devk15 = -9.78 49 REAL(wp) :: devk16 = -48.76 50 REAL(wp) :: devk17 = -14.51 51 REAL(wp) :: devk18 = -23.12 52 REAL(wp) :: devk19 = -26.57 53 REAL(wp) :: devk110 = -29.48 54 ! 55 REAL(wp) :: devk20 = 0.1271 56 REAL(wp) :: devk21 = -0.0219 57 REAL(wp) :: devk22 = 0.1622 58 REAL(wp) :: devk23 = 0.1119 59 REAL(wp) :: devk24 = 0.0466 60 REAL(wp) :: devk25 = -0.0090 61 REAL(wp) :: devk26 = 0.5304 62 REAL(wp) :: devk27 = 0.1211 63 REAL(wp) :: devk28 = 0.1758 64 REAL(wp) :: devk29 = 0.2020 65 REAL(wp) :: devk210 = 0.1622 66 ! 67 REAL(wp) :: devk30 = 0. 68 REAL(wp) :: devk31 = 0. 69 REAL(wp) :: devk32 = 2.608E-3 70 REAL(wp) :: devk33 = -1.409e-3 71 REAL(wp) :: devk34 = 0.316e-3 72 REAL(wp) :: devk35 = -0.942e-3 73 REAL(wp) :: devk36 = 0. 74 REAL(wp) :: devk37 = -0.321e-3 75 REAL(wp) :: devk38 = -2.647e-3 76 REAL(wp) :: devk39 = -3.042e-3 77 REAL(wp) :: devk310 = -2.6080e-3 78 ! 79 REAL(wp) :: devk40 = -3.08E-3 80 REAL(wp) :: devk41 = 1.13E-3 81 REAL(wp) :: devk42 = -2.84E-3 82 REAL(wp) :: devk43 = -5.13E-3 83 REAL(wp) :: devk44 = -4.53e-3 84 REAL(wp) :: devk45 = -3.91e-3 85 REAL(wp) :: devk46 = -11.76e-3 86 REAL(wp) :: devk47 = -2.67e-3 87 REAL(wp) :: devk48 = -5.15e-3 88 REAL(wp) :: devk49 = -4.08e-3 89 REAL(wp) :: devk410 = -2.84e-3 90 ! 91 REAL(wp) :: devk50 = 0.0877E-3 92 REAL(wp) :: devk51 = -0.1475E-3 93 REAL(wp) :: devk52 = 0. 94 REAL(wp) :: devk53 = 0.0794E-3 95 REAL(wp) :: devk54 = 0.09e-3 96 REAL(wp) :: devk55 = 0.054e-3 97 REAL(wp) :: devk56 = 0.3692E-3 98 REAL(wp) :: devk57 = 0.0427e-3 99 REAL(wp) :: devk58 = 0.09e-3 100 REAL(wp) :: devk59 = 0.0714e-3 101 REAL(wp) :: devk510 = 0.0 164 102 165 103 !! $Id$ … … 179 117 INTEGER, INTENT(in) :: kt ! time step 180 118 181 #if ! defined key_sed_off182 119 INTEGER :: ji, jj, ikt 183 REAL(wp) :: ztkel, ztc, ztc2, zpres, ztr 184 REAL(wp) :: zsal, zsal2, zsqrt, zsal15 185 REAL(wp) :: zis, zis2, zisqrt 120 REAL(wp) :: ztc, ztc2 121 REAL(wp) :: zsal, zsal15 186 122 REAL(wp) :: zdens0, zaw, zbw, zcw 187 REAL(wp) :: zbuf1, zbuf2 188 REAL(wp) :: zcpexp, zcpexp2 189 REAL(wp) :: zck1p, zck2p, zck3p, zcksi 190 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zchem_data 191 192 #endif 123 REAL(wp), DIMENSION(jpi,jpj,15) :: zchem_data 193 124 !!---------------------------------------------------------------------- 194 125 195 IF( MOD( kt - 1, nfreq ) /= 0 ) RETURN 196 197 WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt 198 WRITE(numsed,*) ' ' 199 200 201 #if defined key_sed_off 202 CALL sed_chem_off 203 #else 126 127 IF( ln_timing ) CALL timing_start('sed_chem') 128 129 IF (lwp) WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt 130 IF (lwp) WRITE(numsed,*) ' ' 131 204 132 ! reading variables 205 ALLOCATE( zchem_data(jpi,jpj,6) ) ; zchem_data(:,:,:) = 0. 206 207 DO jj = 1,jpj 208 DO ji = 1, jpi 209 ikt = mbkt(ji,jj) 210 IF ( tmask(ji,jj,ikt) == 1 ) THEN 211 zchem_data(ji,jj,1) = ak13 (ji,jj,ikt) 212 zchem_data(ji,jj,2) = ak23 (ji,jj,ikt) 213 zchem_data(ji,jj,3) = akb3 (ji,jj,ikt) 214 zchem_data(ji,jj,4) = akw3 (ji,jj,ikt) 215 zchem_data(ji,jj,5) = aksp (ji,jj,ikt) 216 zchem_data(ji,jj,6) = borat(ji,jj,ikt) 217 ENDIF 133 zchem_data(:,:,:) = rtrn 134 135 IF (ln_sediment_offline) THEN 136 CALL sed_chem_cst 137 ELSE 138 DO jj = 1,jpj 139 DO ji = 1, jpi 140 ikt = mbkt(ji,jj) 141 IF ( tmask(ji,jj,ikt) == 1 ) THEN 142 zchem_data(ji,jj,1) = ak13 (ji,jj,ikt) 143 zchem_data(ji,jj,2) = ak23 (ji,jj,ikt) 144 zchem_data(ji,jj,3) = akb3 (ji,jj,ikt) 145 zchem_data(ji,jj,4) = akw3 (ji,jj,ikt) 146 zchem_data(ji,jj,5) = aksp (ji,jj,ikt) 147 zchem_data(ji,jj,6) = borat (ji,jj,ikt) 148 zchem_data(ji,jj,7) = ak1p3 (ji,jj,ikt) 149 zchem_data(ji,jj,8) = ak2p3 (ji,jj,ikt) 150 zchem_data(ji,jj,9) = ak3p3 (ji,jj,ikt) 151 zchem_data(ji,jj,10)= aksi3 (ji,jj,ikt) 152 zchem_data(ji,jj,11)= sio3eq(ji,jj,ikt) 153 zchem_data(ji,jj,12)= aks3 (ji,jj,ikt) 154 zchem_data(ji,jj,13)= akf3 (ji,jj,ikt) 155 zchem_data(ji,jj,14)= sulfat(ji,jj,ikt) 156 zchem_data(ji,jj,15)= fluorid(ji,jj,ikt) 157 ENDIF 158 ENDDO 218 159 ENDDO 219 ENDDO 220 221 CALL pack_arr ( jpoce, ak1s (1:jpoce), zchem_data(1:jpi,1:jpj,1), iarroce(1:jpoce) ) 222 CALL pack_arr ( jpoce, ak2s (1:jpoce), zchem_data(1:jpi,1:jpj,2), iarroce(1:jpoce) ) 223 CALL pack_arr ( jpoce, akbs (1:jpoce), zchem_data(1:jpi,1:jpj,3), iarroce(1:jpoce) ) 224 CALL pack_arr ( jpoce, akws (1:jpoce), zchem_data(1:jpi,1:jpj,4), iarroce(1:jpoce) ) 225 CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5), iarroce(1:jpoce) ) 226 CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6), iarroce(1:jpoce) ) 227 228 DEALLOCATE( zchem_data ) 160 161 CALL pack_arr ( jpoce, ak1s (1:jpoce), zchem_data(1:jpi,1:jpj,1) , iarroce(1:jpoce) ) 162 CALL pack_arr ( jpoce, ak2s (1:jpoce), zchem_data(1:jpi,1:jpj,2) , iarroce(1:jpoce) ) 163 CALL pack_arr ( jpoce, akbs (1:jpoce), zchem_data(1:jpi,1:jpj,3) , iarroce(1:jpoce) ) 164 CALL pack_arr ( jpoce, akws (1:jpoce), zchem_data(1:jpi,1:jpj,4) , iarroce(1:jpoce) ) 165 CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5) , iarroce(1:jpoce) ) 166 CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6) , iarroce(1:jpoce) ) 167 CALL pack_arr ( jpoce, ak1ps (1:jpoce), zchem_data(1:jpi,1:jpj,7) , iarroce(1:jpoce) ) 168 CALL pack_arr ( jpoce, ak2ps (1:jpoce), zchem_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) ) 169 CALL pack_arr ( jpoce, ak3ps (1:jpoce), zchem_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) ) 170 CALL pack_arr ( jpoce, aksis (1:jpoce), zchem_data(1:jpi,1:jpj,10), iarroce(1:jpoce) ) 171 CALL pack_arr ( jpoce, sieqs (1:jpoce), zchem_data(1:jpi,1:jpj,11), iarroce(1:jpoce) ) 172 CALL pack_arr ( jpoce, aks3s (1:jpoce), zchem_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 173 CALL pack_arr ( jpoce, akf3s (1:jpoce), zchem_data(1:jpi,1:jpj,13), iarroce(1:jpoce) ) 174 CALL pack_arr ( jpoce, sulfats(1:jpoce), zchem_data(1:jpi,1:jpj,14), iarroce(1:jpoce) ) 175 CALL pack_arr ( jpoce, fluorids(1:jpoce), zchem_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 176 ENDIF 229 177 230 178 DO ji = 1, jpoce 231 ztkel = temp(ji) + 273.16232 179 ztc = temp(ji) 233 180 ztc2 = ztc * ztc 234 zpres = press(ji)235 181 ! zqtt = ztkel * 0.01 236 182 zsal = salt(ji) 237 zsal2 = zsal * zsal 238 zsqrt = SQRT( zsal ) 239 zsal15 = zsqrt * zsal 240 zlogt = LOG( ztkel ) 241 ztr = 1./ ztkel 242 ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet) 243 zis = 19.924 * zsal / ( 1000. - 1.005 * zsal ) 244 zis2 = zis * zis 245 zisqrt = SQRT( zis ) 183 zsal15 = SQRT( zsal ) * zsal 246 184 247 185 ! Density of Sea Water - F(temp,sal) [kg/m3] … … 256 194 densSW(ji) = densSW(ji) * 1E-3 ! to get dens in [kg/l] 257 195 258 ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970)259 ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982)260 ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN261 ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.)262 ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP263 ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286264 ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)265 !-----------------------------------------------------------------------------------------266 zcpexp = zpres / ( rgas*ztkel )267 zcpexp2 = zpres * zcpexp268 269 ! For Phodphates (phosphoric acid) (DOE 1994)270 !----------------------------------------------271 zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &272 & + ( cp15*ztr + cp16 ) * zsal273 zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &274 & + ( cp25*ztr + cp26 ) * zsal275 zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) * zsqrt &276 & + ( cp34*ztr + cp35 ) * zsal277 278 ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)279 !--------------------------------------------------------280 zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &281 & + ( cs15*ztr + cs16 ) * zis &282 & + ( cs17*ztr + cs18 ) * zis2 &283 & + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )284 285 286 !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)287 !---------------------------------------------------288 zak1p = EXP ( zck1p )289 zak2p = EXP ( zck2p )290 zak3p = EXP ( zck3p )291 zaksil = EXP ( zcksi )292 293 zbuf1 = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )294 zbuf2 = 0.5 * ( devk4(3) + devk5(3)*ztc )295 aksis(ji) = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )296 297 zbuf1 = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 )298 zbuf2 = 0.5 * ( devk4(6) + devk5(6)*ztc )299 ak1ps(ji) = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )300 301 zbuf1 = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 )302 zbuf2 = 0.5 * ( devk4(7) + devk5(7)*ztc )303 ak2ps(ji) = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )304 305 zbuf1 = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 )306 zbuf2 = 0.5 * ( devk4(8) + devk5(8)*ztc )307 ak3ps(ji) = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )308 309 196 ak12s (ji) = ak1s (ji) * ak2s (ji) 310 197 ak12ps (ji) = ak1ps(ji) * ak2ps(ji) … … 313 200 calcon2(ji) = 0.01028 * ( salt(ji) / 35. ) * densSW(ji) 314 201 ENDDO 315 316 202 317 #endif 203 IF( ln_timing ) CALL timing_stop('sed_chem') 318 204 319 205 END SUBROUTINE sed_chem 320 206 321 #if defined key_sed_off 322 323 SUBROUTINE sed_chem_off 324 !!---------------------------------------------------------------------- 325 !! *** ROUTINE sed_chem_off *** 326 !! 327 !! ** Purpose : compute chemical constants 207 SUBROUTINE ahini_for_at_sed(p_hini) 208 !!--------------------------------------------------------------------- 209 !! *** ROUTINE ahini_for_at *** 328 210 !! 329 !! History : 330 !! ! 04-10 (N. Emprin, M. Gehlen ) Original code 331 !! ! 06-04 (C. Ethe) Re-organization 332 !!---------------------------------------------------------------------- 333 !! * Local declarations 334 INTEGER :: ji 335 336 REAL(wp) :: ztkel, ztc, ztc2, zpres, ztr 337 REAL(wp) :: zsal, zsal2, zsqrt, zsal15 338 REAL(wp) :: zis, zis2, zisqrt 339 REAL(wp) :: zdens0, zaw, zbw, zcw 340 REAL(wp) :: zchl, zst, zft, zbuf1, zbuf2 341 REAL(wp) :: zcpexp, zcpexp2 342 REAL(wp) :: zckb, zck1, zck2, zckw 343 REAL(wp) :: zck1p, zck2p, zck3p, zcksi 344 REAL(wp) :: zak1, zak2, zakb, zakw 345 REAL(wp) :: zaksp0, zksp, zks, zkf 346 211 !! Subroutine returns the root for the 2nd order approximation of the 212 !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 213 !! polynomial) around the local minimum, if it exists. 214 !! Returns * 1E-03_wp if p_alkcb <= 0 215 !! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 216 !! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 217 !! and the 2nd order approximation does not have 218 !! a solution 219 !!--------------------------------------------------------------------- 220 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_hini 221 INTEGER :: ji, jk 222 REAL(wp) :: zca1, zba1 223 REAL(wp) :: zd, zsqrtd, zhmin 224 REAL(wp) :: za2, za1, za0 225 REAL(wp) :: p_dictot, p_bortot, p_alkcb 226 227 IF( ln_timing ) CALL timing_start('ahini_for_at_sed') 228 ! 229 DO jk = 1, jpksed 230 DO ji = 1, jpoce 231 p_alkcb = pwcp(ji,jk,jwalk) / densSW(ji) 232 p_dictot = pwcp(ji,jk,jwdic) / densSW(ji) 233 p_bortot = borats(ji) / densSW(ji) 234 IF (p_alkcb <= 0.) THEN 235 p_hini(ji,jk) = 1.e-3 236 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 237 p_hini(ji,jk) = 1.e-10_wp 238 ELSE 239 zca1 = p_dictot/( p_alkcb + rtrn ) 240 zba1 = p_bortot/ (p_alkcb + rtrn ) 241 ! Coefficients of the cubic polynomial 242 za2 = aKbs(ji)*(1. - zba1) + ak1s(ji)*(1.-zca1) 243 za1 = ak1s(ji)*akbs(ji)*(1. - zba1 - zca1) & 244 & + ak1s(ji)*ak2s(ji)*(1. - (zca1+zca1)) 245 za0 = ak1s(ji)*ak2s(ji)*akbs(ji)*(1. - zba1 - (zca1+zca1)) 246 ! Taylor expansion around the minimum 247 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 248 ! for the minimum close to the root 249 250 IF(zd > 0.) THEN ! If the discriminant is positive 251 zsqrtd = SQRT(zd) 252 IF(za2 < 0) THEN 253 zhmin = (-za2 + zsqrtd)/3. 254 ELSE 255 zhmin = -za1/(za2 + zsqrtd) 256 ENDIF 257 p_hini(ji,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 258 ELSE 259 p_hini(ji,jk) = 1.e-7 260 ENDIF 261 ! 262 ENDIF 263 END DO 264 END DO 265 ! 266 IF( ln_timing ) CALL timing_stop('ahini_for_at_sed') 267 ! 268 END SUBROUTINE ahini_for_at_sed 269 270 !=============================================================================== 271 SUBROUTINE anw_infsup_sed( p_alknw_inf, p_alknw_sup ) 272 273 ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 274 ! contributions to total alkalinity (the infimum and the supremum), i.e 275 ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 276 277 ! Argument variables 278 INTEGER :: jk 279 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_inf 280 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_sup 281 282 DO jk = 1, jpksed 283 p_alknw_inf(:,jk) = -pwcp(:,jk,jwpo4) / densSW(:) 284 p_alknw_sup(:,jk) = (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil) & 285 & + borats(:) ) / densSW(:) 286 END DO 287 288 END SUBROUTINE anw_infsup_sed 289 290 291 SUBROUTINE solve_at_general_sed( p_hini, zhi ) 292 293 ! Universal pH solver that converges from any given initial value, 294 ! determines upper an lower bounds for the solution if required 295 296 ! Argument variables 297 !-------------------- 298 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(IN) :: p_hini 299 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: zhi 300 301 ! Local variables 302 !----------------- 303 INTEGER :: ji, jk, jn 304 REAL(wp) :: zh_ini, zh, zh_prev, zh_lnfactor 305 REAL(wp) :: zdelta, zh_delta 306 REAL(wp) :: zeqn, zdeqndh, zalka 307 REAL(wp) :: aphscale 308 REAL(wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 309 REAL(wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 310 REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 311 REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 312 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 313 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 314 REAL(wp) :: zalk_wat, zdalk_wat 315 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 316 LOGICAL :: l_exitnow 317 REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 318 REAL(wp), DIMENSION(jpoce,jpksed) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 319 320 IF( ln_timing ) CALL timing_start('solve_at_general_sed') 321 ! Allocate temporary workspace 322 CALL anw_infsup_sed( zalknw_inf, zalknw_sup ) 323 324 rmask(:,:) = 1.0 325 zhi(:,:) = 0. 326 327 ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 328 DO jk = 1, jpksed 329 DO ji = 1, jpoce 330 IF (rmask(ji,jk) == 1.) THEN 331 p_alktot = pwcp(ji,jk,jwalk) / densSW(ji) 332 aphscale = 1. + sulfats(ji)/aks3s(ji) 333 zh_ini = p_hini(ji,jk) 334 335 zdelta = (p_alktot-zalknw_inf(ji,jk))**2 + 4.*akws(ji) / aphscale 336 337 IF(p_alktot >= zalknw_inf(ji,jk)) THEN 338 zh_min(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_inf(ji,jk) + SQRT(zdelta) ) 339 ELSE 340 zh_min(ji,jk) = aphscale * (-(p_alktot-zalknw_inf(ji,jk)) + SQRT(zdelta) ) / 2. 341 ENDIF 342 343 zdelta = (p_alktot-zalknw_sup(ji,jk))**2 + 4.*akws(ji) / aphscale 344 345 IF(p_alktot <= zalknw_sup(ji,jk)) THEN 346 zh_max(ji,jk) = aphscale * (-(p_alktot-zalknw_sup(ji,jk)) + SQRT(zdelta) ) / 2. 347 ELSE 348 zh_max(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_sup(ji,jk) + SQRT(zdelta) ) 349 ENDIF 350 351 zhi(ji,jk) = MAX(MIN(zh_max(ji,jk), zh_ini), zh_min(ji,jk)) 352 ENDIF 353 END DO 354 END DO 355 356 zeqn_absmin(:,:) = HUGE(1._wp) 357 358 DO jn = 1, jp_maxniter_atgen 359 DO jk = 1, jpksed 360 DO ji = 1, jpoce 361 IF (rmask(ji,jk) == 1.) THEN 362 363 p_alktot = pwcp(ji,jk,jwalk) / densSW(ji) 364 zdic = pwcp(ji,jk,jwdic) / densSW(ji) 365 zbot = borats(ji) / densSW(ji) 366 zpt = pwcp(ji,jk,jwpo4) / densSW(ji) 367 zsit = pwcp(ji,jk,jwsil) / densSW(ji) 368 zst = sulfats(ji) 369 zft = fluorids(ji) 370 aphscale = 1. + sulfats(ji)/aks3s(ji) 371 zh = zhi(ji,jk) 372 zh_prev = zh 373 374 ! H2CO3 - HCO3 - CO3 : n=2, m=0 375 znumer_dic = 2.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji) 376 zdenom_dic = ak1s(ji)*ak2s(ji) + zh*(ak1s(ji) + zh) 377 zalk_dic = zdic * (znumer_dic/zdenom_dic) 378 zdnumer_dic = ak1s(ji)*ak1s(ji)*ak2s(ji) + zh & 379 *(4.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji)) 380 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 381 382 383 ! B(OH)3 - B(OH)4 : n=1, m=0 384 znumer_bor = akbs(ji) 385 zdenom_bor = akbs(ji) + zh 386 zalk_bor = zbot * (znumer_bor/zdenom_bor) 387 zdnumer_bor = akbs(ji) 388 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 389 390 391 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 392 znumer_po4 = 3.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji) & 393 & + zh*(2.*ak1ps(ji)*ak2ps(ji) + zh* ak1ps(ji)) 394 zdenom_po4 = ak1ps(ji)*ak2ps(ji)*ak3ps(ji) & 395 & + zh*( ak1ps(ji)*ak2ps(ji) + zh*(ak1ps(ji) + zh)) 396 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 397 zdnumer_po4 = ak1ps(ji)*ak2ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji) & 398 & + zh*(4.*ak1ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji) & 399 & + zh*(9.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji) & 400 & + ak1ps(ji)*ak1ps(ji)*ak2ps(ji) & 401 & + zh*(4.*ak1ps(ji)*ak2ps(ji) + zh * ak1ps(ji) ) ) ) 402 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 403 404 ! H4SiO4 - H3SiO4 : n=1, m=0 405 znumer_sil = aksis(ji) 406 zdenom_sil = aksis(ji) + zh 407 zalk_sil = zsit * (znumer_sil/zdenom_sil) 408 zdnumer_sil = aksis(ji) 409 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 410 411 ! HSO4 - SO4 : n=1, m=1 412 aphscale = 1.0 + zst/aks3s(ji) 413 znumer_so4 = aks3s(ji) * aphscale 414 zdenom_so4 = aks3s(ji) * aphscale + zh 415 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 416 zdnumer_so4 = aks3s(ji) * aphscale 417 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 418 419 ! HF - F : n=1, m=1 420 znumer_flu = akf3s(ji) 421 zdenom_flu = akf3s(ji) + zh 422 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 423 zdnumer_flu = akf3s(ji) 424 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 425 426 ! H2O - OH 427 zalk_wat = akws(ji)/zh - zh/aphscale 428 zdalk_wat = -akws(ji)/zh**2 - 1./aphscale 429 430 ! CALCULATE [ALK]([CO3--], [HCO3-]) 431 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 432 & + zalk_so4 + zalk_flu & 433 & + zalk_wat - p_alktot 434 435 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 436 & + zalk_so4 + zalk_flu + zalk_wat) 437 438 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 439 & + zdalk_so4 + zdalk_flu + zdalk_wat 440 441 ! Adapt bracketing interval 442 IF(zeqn > 0._wp) THEN 443 zh_min(ji,jk) = zh_prev 444 ELSEIF(zeqn < 0._wp) THEN 445 zh_max(ji,jk) = zh_prev 446 ENDIF 447 448 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jk)) THEN 449 ! if the function evaluation at the current point is 450 ! not decreasing faster than with a bisection step (at least linearly) 451 ! in absolute value take one bisection step on [ph_min, ph_max] 452 ! ph_new = (ph_min + ph_max)/2d0 453 ! 454 ! In terms of [H]_new: 455 ! [H]_new = 10**(-ph_new) 456 ! = 10**(-(ph_min + ph_max)/2d0) 457 ! = SQRT(10**(-(ph_min + phmax))) 458 ! = SQRT(zh_max * zh_min) 459 zh = SQRT(zh_max(ji,jk) * zh_min(ji,jk)) 460 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 461 ELSE 462 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 463 ! = -zdeqndh * LOG(10) * [H] 464 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 465 ! 466 ! pH_new = pH_old + \deltapH 467 ! 468 ! [H]_new = 10**(-pH_new) 469 ! = 10**(-pH_old - \Delta pH) 470 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 471 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 472 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 473 474 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 475 476 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 477 zh = zh_prev*EXP(zh_lnfactor) 478 ELSE 479 zh_delta = zh_lnfactor*zh_prev 480 zh = zh_prev + zh_delta 481 ENDIF 482 483 IF( zh < zh_min(ji,jk) ) THEN 484 ! if [H]_new < [H]_min 485 ! i.e., if ph_new > ph_max then 486 ! take one bisection step on [ph_prev, ph_max] 487 ! ph_new = (ph_prev + ph_max)/2d0 488 ! In terms of [H]_new: 489 ! [H]_new = 10**(-ph_new) 490 ! = 10**(-(ph_prev + ph_max)/2d0) 491 ! = SQRT(10**(-(ph_prev + phmax))) 492 ! = SQRT([H]_old*10**(-ph_max)) 493 ! = SQRT([H]_old * zh_min) 494 zh = SQRT(zh_prev * zh_min(ji,jk)) 495 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 496 ENDIF 497 498 IF( zh > zh_max(ji,jk) ) THEN 499 ! if [H]_new > [H]_max 500 ! i.e., if ph_new < ph_min, then 501 ! take one bisection step on [ph_min, ph_prev] 502 ! ph_new = (ph_prev + ph_min)/2d0 503 ! In terms of [H]_new: 504 ! [H]_new = 10**(-ph_new) 505 ! = 10**(-(ph_prev + ph_min)/2d0) 506 ! = SQRT(10**(-(ph_prev + ph_min))) 507 ! = SQRT([H]_old*10**(-ph_min)) 508 ! = SQRT([H]_old * zhmax) 509 zh = SQRT(zh_prev * zh_max(ji,jk)) 510 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 511 ENDIF 512 ENDIF 513 514 zeqn_absmin(ji,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jk)) 515 516 ! Stop iterations once |\delta{[H]}/[H]| < rdel 517 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 518 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 519 ! Alternatively: 520 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 521 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 522 ! < 1/LOG(10) * rdel 523 524 ! Hence |zeqn/(zdeqndh*zh)| < rdel 525 526 ! rdel <-- pp_rdel_ah_target 527 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 528 529 IF(l_exitnow) THEN 530 rmask(ji,jk) = 0. 531 ENDIF 532 533 zhi(ji,jk) = zh 534 535 IF(jn >= jp_maxniter_atgen) THEN 536 zhi(ji,jk) = -1._wp 537 ENDIF 538 539 ENDIF 540 END DO 541 END DO 542 END DO 543 ! 544 IF( ln_timing ) CALL timing_stop('solve_at_general_sed') 545 546 END SUBROUTINE solve_at_general_sed 547 548 SUBROUTINE sed_chem_cst 549 !!--------------------------------------------------------------------- 550 !! *** ROUTINE sed_chem_cst *** 551 !! 552 !! ** Purpose : Sea water chemistry computed following MOCSY protocol 553 !! Computation is done at the bottom of the ocean only 554 !! 555 !! ** Method : - ... 556 !!--------------------------------------------------------------------- 557 INTEGER :: ji 558 REAL(wp), DIMENSION(jpoce) :: saltprac, temps 559 REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2 560 REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 561 REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2 562 REAL(wp) :: zsqrt, ztr , zlogt , zcek1, zc1, zplat 563 REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1, za2 564 REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw 565 REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 566 REAL(wp) :: zst , zft , zcks , zckf , zaksp1 567 REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total 568 !!--------------------------------------------------------------------- 569 ! 570 IF( ln_timing ) CALL timing_start('sed_chem_cst') 571 ! 572 ! Computation of chemical constants require practical salinity 573 ! Thus, when TEOS08 is used, absolute salinity is converted to 574 ! practical salinity 575 ! ------------------------------------------------------------- 576 IF (neos == -1) THEN 577 saltprac(:) = salt(:) * 35.0 / 35.16504 578 ELSE 579 saltprac(:) = temp(:) 580 ENDIF 581 582 ! 583 ! Computations of chemical constants require in situ temperature 584 ! Here a quite simple formulation is used to convert 585 ! potential temperature to in situ temperature. The errors is less than 586 ! 0.04°C relative to an exact computation 587 ! --------------------------------------------------------------------- 588 DO ji = 1, jpoce 589 zpres = zkbot(ji) / 1000. 590 za1 = 0.04 * ( 1.0 + 0.185 * temp(ji) + 0.035 * (saltprac(ji) - 35.0) ) 591 za2 = 0.0075 * ( 1.0 - temp(ji) / 30.0 ) 592 temps(ji) = temp(ji) - za1 * zpres + za2 * zpres**2 593 END DO 347 594 348 595 ! CHEMICAL CONSTANTS - DEEP OCEAN 349 !------------------------------------- 350 ! [chem constants]=mol/kg solution (or (mol/kg sol)2 for akws and aksp) 351 596 ! ------------------------------- 352 597 DO ji = 1, jpoce 353 ztkel = temp(ji) + 273.16 354 ztc = temp(ji) 355 ztc2 = ztc * ztc 356 zpres = press(ji) 357 ! zqtt = ztkel * 0.01 358 zsal = salt(ji) 359 zsal2 = zsal * zsal 360 zsqrt = SQRT( zsal ) 598 ! SET PRESSION ACCORDING TO SAUNDER (1980) 599 zc1 = 5.92E-3 600 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*zkbot(ji)))) / 4.42E-6 601 zpres = zpres / 10.0 602 603 ! SET ABSOLUTE TEMPERATURE 604 ztkel = temps(ji) + 273.15 605 zsal = saltprac(ji) 606 zsqrt = SQRT( zsal ) 361 607 zsal15 = zsqrt * zsal 362 zlogt = LOG( ztkel ) 363 ztr = 1./ ztkel 364 ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet) 365 zis = 19.924 * zsal / ( 1000. - 1.005 * zsal ) 366 zis2 = zis * zis 367 zisqrt = SQRT( zis ) 368 369 370 ! Density of Sea Water - F(temp,sal) [kg/m3] 371 zdens0 = Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 & 372 + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 & 373 + Ddsw(6) * ztc * ztc2 * ztc2 374 zaw = Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 & 375 + Adsw(5) * ztc2 * ztc2 376 zbw = Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2 377 zcw = Cdsw 378 densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal 379 densSW(ji) = densSW(ji) * 1E-3 ! to get dens in [kg/l] 380 381 382 ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970) 383 ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982) 384 ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN 385 ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.) 386 ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP 387 ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286 388 ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285) 389 !----------------------------------------------------------------------------------------- 390 zcpexp = zpres / ( rgas*ztkel ) 608 zlogt = LOG( ztkel ) 609 ztr = 1. / ztkel 610 zis = 19.924 * zsal / ( 1000.- 1.005 * zsal ) 611 zis2 = zis * zis 612 zisqrt = SQRT( zis ) 613 ztc = temps(ji) 614 615 ! CHLORINITY (WOOSTER ET AL., 1969) 616 zcl = zsal / 1.80655 617 618 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 619 zst = 0.14 * zcl /96.062 620 621 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 622 zft = 0.000067 * zcl /18.9984 623 624 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 625 zcks = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt & 626 & + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 627 & + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis & 628 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 629 & + LOG(1.0 - 0.001005 * zsal)) 630 631 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 632 zckf = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt & 633 & + LOG(1.0d0 - 0.001005d0*zsal) & 634 & + LOG(1.0d0 + zst/zcks)) 635 636 ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 637 zckb= (-8966.90 - 2890.53*zsqrt - 77.942*zsal & 638 & + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr & 639 & + (148.0248 + 137.1942*zsqrt + 1.62142*zsal) & 640 & + (-24.4344 - 25.085*zsqrt - 0.2474*zsal) & 641 & * zlogt + 0.053105*zsqrt*ztkel 642 643 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 644 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 645 zck1 = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt & 646 - 0.011555*zsal + 0.0001152*zsal*zsal) 647 zck2 = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt & 648 - 0.01781*zsal + 0.0001122*zsal*zsal) 649 650 ! PKW (H2O) (MILLERO, 1995) from composite data 651 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 652 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 653 654 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 655 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 656 & + (-106.736*ztr + 0.69171) * zsqrt & 657 & + (-0.65643*ztr - 0.01844) * zsal 658 659 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 660 & + (-160.340*ztr + 1.3566)*zsqrt & 661 & + (0.37335*ztr - 0.05778)*zsal 662 663 zck3p = -3070.75*ztr - 18.126 & 664 & + (17.27039*ztr + 2.81197) * zsqrt & 665 & + (-44.99486*ztr - 0.09984) * zsal 666 667 ! CONSTANT FOR SILICATE, MILLERO (1995) 668 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 669 & + (-458.79*ztr + 3.5913) * zisqrt & 670 & + (188.74*ztr - 1.5998) * zis & 671 & + (-12.1652*ztr + 0.07871) * zis2 & 672 & + LOG(1.0 - 0.001005*zsal) 673 674 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 675 ! (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 676 zaksp0 = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel ) & 677 & + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt & 678 & - 0.07711*zsal + 0.0041249*zsal15 679 680 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 681 zak1 = 10**(zck1) * total2SWS 682 zak2 = 10**(zck2) * total2SWS 683 zakb = EXP( zckb ) * total2SWS 684 zakw = EXP( zckw ) 685 zaksp1 = 10**(zaksp0) 686 zak1p = exp( zck1p ) 687 zak2p = exp( zck2p ) 688 zak3p = exp( zck3p ) 689 zaksi = exp( zcksi ) 690 zckf = zckf * total2SWS 691 692 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 693 ! (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE 694 ! IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 695 ! TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN 696 ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 697 ! MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 698 ! WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 699 ! & GIESKES (1970), P. 1285-1286 (THE SMALL 700 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 701 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 702 zcpexp = zpres / (rgas*ztkel) 391 703 zcpexp2 = zpres * zcpexp 392 704 393 394 ! chlorinity (WOOSTER ET AL., 1969) 395 !--------------------------------------- 396 zchl = zsal * salchl 397 398 ! total sulfate concentration [mol/kg soln] 399 ! -------------------------------------- 400 zst = st1 * zchl * st2 401 402 ! total fluoride concentration [mol/kg soln] 403 ! -------------------------------------- 404 zft = ft1 * zchl * ft2 405 406 ! dissociation constant for carbonate (Mehrback 74 - Dickson & Millero 87) 407 !--------------------------------------------------------------------------- 408 zck1 = c10*ztr - c11 + c12*zlogt - c13*zsal + c14*zsal2 409 zck2 = c20*ztr + c21 - c22*zlogt - c23*zsal + c24*zsal2 410 411 ! dissociation constant for sulfates (Dickson 1990) 412 !-------------------------------------------------- 413 zks = EXP( ks0 + ks1*ztr + ks2*zlogt & 414 & + ( ks3*ztr + ks4 + ks5*zlogt ) * zisqrt & 415 & + ( ks6*ztr + ks7 + ks8*zlogt ) * zis & 416 & + ks9*ztr*zis*zisqrt + ks10*ztr*zis2 & 417 & + LOG( ks11 + ks12*zsal ) ) 418 419 ! dissociation constant for fluorides (Dickson and Riley 79) 420 !-------------------------------------------------- 421 zkf = EXP( kf0 + kf1*ztr + kf2*zisqrt + LOG( kf3 + kf4*zsal ) ) 422 423 ! dissociation constant for borates (Doe 94) 424 !-------------------------------------------------- 425 zckb = ( cb0 + cb1*zsqrt + cb2*zsal + cb3*zsal15 + cb4*zsal2) * ztr & 426 & + ( cb5 + cb6*zsqrt + cb7*zsal) & 427 & + ( cb8 + cb9*zsqrt + cb10*zsal) * zlogt & 428 & + cb11*zsqrt*ztkel + LOG( ( 1. + zst/zks + zft/zkf ) / ( 1. + zst/zks ) ) 429 430 ! PKW (H2O) (DICKSON AND RILEY, 1979) 431 !-------------------------------------- 432 zckw = cw0*ztr + cw1 + cw2*zlogt & 433 & +( cw3*ztr + cw4 + cw5*zlogt )* zsqrt + cw6*zsal 434 435 ! For Phodphates (phosphoric acid) (DOE 1994) 436 !---------------------------------------------- 437 zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt & 438 & + ( cp15*ztr + cp16 ) * zsal 439 zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt & 440 & + ( cp25*ztr + cp26 ) * zsal 441 zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) * zsqrt & 442 & + ( cp34*ztr + cp35 ) * zsal 443 444 ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP) 445 !-------------------------------------------------------- 446 zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt & 447 & + ( cs15*ztr + cs16 ) * zis & 448 & + ( cs17*ztr + cs18 ) * zis2 & 449 & + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal ) 450 451 ! apparent solublity product K'SP of calcite in seawater 452 ! (S=27-43, T=2-25 DEG C) AT pres =0 (INGLE, 1975, EQ. 6) 453 ! prob: olivier a log = ln et C. Heize a LOG10(sal) 454 ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)+akcc4*tkel*tkel) 455 ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)+akcc4*tkel*tkel) 456 !-------------------------------------------------------------------- 457 zaksp0 = akcc1 + akcc2*ztkel + akcc3*ztr + akcc4 * LOG10(ztkel) & 458 & + ( akcc5 + akcc6*ztkel+ akcc7*ztr ) * zsqrt & 459 & + akcc8*zsal + akcc9*zsal15 460 461 !K1, K2 of carbonic acid, KB of boric acid, KW (H2O) 462 !--------------------------------------------------- 463 zak1 = 10**( -zck1 ) 464 zak2 = 10**( -zck2 ) 465 zakb = EXP ( zckb ) 466 zakw = EXP ( zckw ) 467 zksp = 10**( zaksp0 ) 468 469 470 471 ! KB of boric acid, K1,K2 of carbonic acid pressure correction 472 ! after Culberson and AND Pytkowicz (1968) (CF. BROECKER ET AL., 1982) Millero 95 473 !-------------------------------------------------------------------------------- 474 zbuf1 = - ( devk1(1) + devk2(1)*ztc + devk3(1)*ztc2 ) 475 zbuf2 = 0.5 * ( devk4(1) + devk5(1)*ztc ) 476 ak1s(ji) = zak1 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 477 478 zbuf1 = -( devk1(2) + devk2(2)*ztc + devk3(2)*ztc2 ) 479 zbuf2 = 0.5 * ( devk4(2) + devk5(2)*ztc ) 480 ak2s(ji) = zak2 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 481 482 zbuf1 = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 ) 483 zbuf2 = 0.5 * ( devk4(3) + devk5(3) * ztc ) 484 akbs(ji) = zakb * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 485 486 zbuf1 = - ( devk1(4) + devk2(4)*ztc + devk3(4)*ztc2 ) 487 zbuf2 = 0.5 * ( devk4(4) + devk5(4)*ztc ) 488 akws(ji) = zakw * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 489 490 491 ! APPARENT SOLUBILITY PRODUCT K''SP OF CALCITE (OR ARAGONITE) 492 ! AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970) 493 ! (P. 1285) AND BERNER (1976) 494 !----------------------------------------------------------------- 495 ! aksp(ji) = aksp0*exp(zcpexp*(devks-devkst*tc)) 496 ! or following Mucci 497 zbuf1 = - ( devk1(5) + devk2(5)*ztc + devk3(5)*ztc2 ) 498 zbuf2 = 0.5 *( devk4(5) + devk5(5)*ztc ) 499 aksps(ji) = zksp * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 500 501 ! For Phodphates (phosphoric acid) (DOE 1994) 502 !---------------------------------------------- 503 zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt & 504 & + ( cp15*ztr + cp16 ) * zsal 505 zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt & 506 & + ( cp25*ztr + cp26 ) * zsal 507 zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) * zsqrt & 508 & + ( cp34*ztr + cp35 ) * zsal 509 510 ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP) 511 !-------------------------------------------------------- 512 zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt & 513 & + ( cs15*ztr + cs16 ) * zis & 514 & + ( cs17*ztr + cs18 ) * zis2 & 515 & + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal ) 516 517 518 !K1, K2 of carbonic acid, KB of boric acid, KW (H2O) 519 !--------------------------------------------------- 520 zak1p = EXP ( zck1p ) 521 zak2p = EXP ( zck2p ) 522 zak3p = EXP ( zck3p ) 523 zaksil = EXP ( zcksi ) 524 525 zbuf1 = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 ) 526 zbuf2 = 0.5 * ( devk4(3) + devk5(3)*ztc ) 527 aksis(ji) = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 528 529 zbuf1 = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 ) 530 zbuf2 = 0.5 * ( devk4(6) + devk5(6)*ztc ) 531 ak1ps(ji) = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 532 533 zbuf1 = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 ) 534 zbuf2 = 0.5 * ( devk4(7) + devk5(7)*ztc ) 535 ak2ps(ji) = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 536 537 zbuf1 = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 ) 538 zbuf2 = 0.5 * ( devk4(8) + devk5(8)*ztc ) 539 ak3ps(ji) = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 ) 540 541 ! total borat concentration. [mol/l] 542 ! or from Millero 1995 [mol/l] : borat(l) = 0.000416_8*(sal/35._8)*densSW(l) 543 ! -------------------------------------------------------------------------- 544 borats(ji) = bor1 * zchl * bor2 * densSW(ji) 545 546 ak12s (ji) = ak1s (ji) * ak2s (ji) 547 ak12ps (ji) = ak1ps(ji) * ak2ps(ji) 548 ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji) 549 550 calcon2(ji) = 0.01028 * ( zsal / 35. ) * densSW(ji) 551 552 ENDDO 553 554 END SUBROUTINE sed_chem_off 555 556 #endif 557 558 #else 559 !!====================================================================== 560 !! MODULE sedchem : Dummy module 561 !!====================================================================== 562 !! $Id$ 563 CONTAINS 564 SUBROUTINE sed_chem( kt ) ! Empty routine 565 INTEGER, INTENT(in) :: kt 566 WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt 567 END SUBROUTINE sed_chem 568 569 !!====================================================================== 570 571 #endif 705 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 706 ! CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 707 ! (CF. BROECKER ET AL., 1982) 708 709 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 710 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 711 ak1s(ji) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 712 713 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 714 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 715 ak2s(ji) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 716 717 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 718 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 719 akbs(ji) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 720 721 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 722 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 723 akws(ji) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 724 725 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 726 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 727 aks3s(ji) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 728 729 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 730 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 731 akf3s(ji) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 732 733 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 734 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 735 ak1ps(ji) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 736 737 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 738 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 739 ak2ps(ji) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 740 741 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 742 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 743 aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 744 745 ! CONVERT FROM DIFFERENT PH SCALES 746 total2free = 1.0/(1.0 + zst/aks3s(ji)) 747 free2SWS = 1. + zst/aks3s(ji) + zft/akf3s(ji) 748 total2SWS = total2free * free2SWS 749 SWS2total = 1.0 / total2SWS 750 751 ! Convert to total scale 752 ak1s(ji) = ak1s(ji) * SWS2total 753 ak2s(ji) = ak2s(ji) * SWS2total 754 akbs(ji) = akbs(ji) * SWS2total 755 akws(ji) = akws(ji) * SWS2total 756 ak1ps(ji) = ak1ps(ji) * SWS2total 757 ak2ps(ji) = ak2ps(ji) * SWS2total 758 ak3ps(ji) = ak3ps(ji) * SWS2total 759 aksis(ji) = aksis(ji) * SWS2total 760 akf3s(ji) = akf3s(ji) / total2free 761 762 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 763 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 764 ! (P. 1285) AND BERNER (1976) 765 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 766 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 767 aksps(ji) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 768 769 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 770 borats(ji) = 0.0002414 * zcl / 10.811 771 sulfats(ji) = zst 772 fluorids(ji) = zft 773 774 ! Iron and SIO3 saturation concentration from ... 775 sieqs(ji) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 776 END DO 777 ! 778 IF( ln_timing ) CALL timing_stop('sed_chem_cst') 779 ! 780 END SUBROUTINE sed_chem_cst 781 572 782 573 783 END MODULE sedchem
Note: See TracChangeset
for help on using the changeset viewer.