Changeset 7607 for branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC
- Timestamp:
- 2017-01-25T16:37:31+01:00 (7 years ago)
- Location:
- branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
r6943 r7607 11 11 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 12 12 !! ! 2011-02 (J. Simeon, J.Orr ) update O2 solubility constants 13 !! 3.6 ! 2016-03 (O. Aumont) Change chemistry to MOCSY standards 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_pisces 15 16 !!---------------------------------------------------------------------- 16 !! 'key_pisces 'PISCES bio-model17 !! 'key_pisces*' PISCES bio-model 17 18 !!---------------------------------------------------------------------- 18 19 !! p4z_che : Sea water chemistry computed following OCMIP protocol … … 22 23 USE sms_pisces ! PISCES Source Minus Sink variables 23 24 USE lib_mpp ! MPP library 25 USE eosbn2, ONLY : nn_eos 24 26 25 27 IMPLICIT NONE 26 28 PRIVATE 27 29 28 PUBLIC p4z_che ! 29 PUBLIC p4z_che_alloc ! 30 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2 30 PUBLIC p4z_che ! 31 PUBLIC p4z_che_alloc ! 32 PUBLIC p4z_che_ahini ! 33 PUBLIC p4z_che_solve_hi ! 34 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol ! solubility of Fe 35 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tempis ! In situ temperature 41 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: salinprac ! Practical salinity 42 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ??? 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ??? 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akf3 !: ??? 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aks3 !: ??? 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak1p3 !: ??? 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak2p3 !: ??? 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak3p3 !: ??? 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksi3 !: ??? 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ??? 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fluorid !: ??? 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sulfat !: ??? 54 55 !!* Variable for chemistry of the CO2 cycle 36 56 37 57 REAL(wp), PUBLIC :: atcox = 0.20946 ! units atm 38 58 39 REAL(wp) :: salchl = 1. / 1.80655 ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)40 59 REAL(wp) :: o2atm = 1. / ( 1000. * 0.20946 ) 41 60 42 REAL(wp) :: rgas = 83.14472 ! universal gas constants 43 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 44 45 REAL(wp) :: bor1 = 0.00023 ! borat constants 46 REAL(wp) :: bor2 = 1. / 10.82 47 48 REAL(wp) :: st1 = 0.14 ! constants for calculate concentrations for sulfate 49 REAL(wp) :: st2 = 1./96.062 ! (Morris & Riley 1966) 50 51 REAL(wp) :: ft1 = 0.000067 ! constants for calculate concentrations for fluorides 52 REAL(wp) :: ft2 = 1./18.9984 ! (Dickson & Riley 1979 ) 53 54 ! ! volumetric solubility constants for o2 in ml/L 55 REAL(wp) :: ox0 = 2.00856 ! from Table 1 for Eq 8 of Garcia and Gordon, 1992. 56 REAL(wp) :: ox1 = 3.22400 ! corrects for moisture and fugacity, but not total atmospheric pressure 57 REAL(wp) :: ox2 = 3.99063 ! Original PISCES code noted this was a solubility, but 58 REAL(wp) :: ox3 = 4.80299 ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2) 59 REAL(wp) :: ox4 = 9.78188e-1 ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2 60 REAL(wp) :: ox5 = 1.71069 ! and atcox = 0.20946 to add the 1/atm dimension. 61 REAL(wp) :: ox6 = -6.24097e-3 62 REAL(wp) :: ox7 = -6.93498e-3 63 REAL(wp) :: ox8 = -6.90358e-3 64 REAL(wp) :: ox9 = -4.29155e-3 65 REAL(wp) :: ox10 = -3.11680e-7 61 REAL(wp) :: rgas = 83.14472 ! universal gas constants 62 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 66 63 67 64 ! ! coeff. for seawater pressure correction : millero 95 68 65 ! ! AGRIF doesn't like the DATA instruction 69 REAL(wp) :: devk11 = -25.5 70 REAL(wp) :: devk12 = -15.82 71 REAL(wp) :: devk13 = -29.48 72 REAL(wp) :: devk14 = -25.60 73 REAL(wp) :: devk15 = -48.76 66 REAL(wp) :: devk10 = -25.5 67 REAL(wp) :: devk11 = -15.82 68 REAL(wp) :: devk12 = -29.48 69 REAL(wp) :: devk13 = -20.02 70 REAL(wp) :: devk14 = -18.03 71 REAL(wp) :: devk15 = -9.78 72 REAL(wp) :: devk16 = -48.76 73 REAL(wp) :: devk17 = -14.51 74 REAL(wp) :: devk18 = -23.12 75 REAL(wp) :: devk19 = -26.57 76 REAL(wp) :: devk110 = -29.48 74 77 ! 75 REAL(wp) :: devk21 = 0.1271 76 REAL(wp) :: devk22 = -0.0219 77 REAL(wp) :: devk23 = 0.1622 78 REAL(wp) :: devk24 = 0.2324 79 REAL(wp) :: devk25 = 0.5304 78 REAL(wp) :: devk20 = 0.1271 79 REAL(wp) :: devk21 = -0.0219 80 REAL(wp) :: devk22 = 0.1622 81 REAL(wp) :: devk23 = 0.1119 82 REAL(wp) :: devk24 = 0.0466 83 REAL(wp) :: devk25 = -0.0090 84 REAL(wp) :: devk26 = 0.5304 85 REAL(wp) :: devk27 = 0.1211 86 REAL(wp) :: devk28 = 0.1758 87 REAL(wp) :: devk29 = 0.2020 88 REAL(wp) :: devk210 = 0.1622 80 89 ! 90 REAL(wp) :: devk30 = 0. 81 91 REAL(wp) :: devk31 = 0. 82 REAL(wp) :: devk32 = 0. 83 REAL(wp) :: devk33 = 2.608E-3 84 REAL(wp) :: devk34 = -3.6246E-3 85 REAL(wp) :: devk35 = 0. 92 REAL(wp) :: devk32 = 2.608E-3 93 REAL(wp) :: devk33 = -1.409e-3 94 REAL(wp) :: devk34 = 0.316e-3 95 REAL(wp) :: devk35 = -0.942e-3 96 REAL(wp) :: devk36 = 0. 97 REAL(wp) :: devk37 = -0.321e-3 98 REAL(wp) :: devk38 = -2.647e-3 99 REAL(wp) :: devk39 = -3.042e-3 100 REAL(wp) :: devk310 = -2.6080e-3 86 101 ! 87 REAL(wp) :: devk41 = -3.08E-3 88 REAL(wp) :: devk42 = 1.13E-3 89 REAL(wp) :: devk43 = -2.84E-3 90 REAL(wp) :: devk44 = -5.13E-3 91 REAL(wp) :: devk45 = -11.76E-3 102 REAL(wp) :: devk40 = -3.08E-3 103 REAL(wp) :: devk41 = 1.13E-3 104 REAL(wp) :: devk42 = -2.84E-3 105 REAL(wp) :: devk43 = -5.13E-3 106 REAL(wp) :: devk44 = -4.53e-3 107 REAL(wp) :: devk45 = -3.91e-3 108 REAL(wp) :: devk46 = -11.76e-3 109 REAL(wp) :: devk47 = -2.67e-3 110 REAL(wp) :: devk48 = -5.15e-3 111 REAL(wp) :: devk49 = -4.08e-3 112 REAL(wp) :: devk410 = -2.84e-3 92 113 ! 93 REAL(wp) :: devk51 = 0.0877E-3 94 REAL(wp) :: devk52 = -0.1475E-3 95 REAL(wp) :: devk53 = 0. 96 REAL(wp) :: devk54 = 0.0794E-3 97 REAL(wp) :: devk55 = 0.3692E-3 114 REAL(wp) :: devk50 = 0.0877E-3 115 REAL(wp) :: devk51 = -0.1475E-3 116 REAL(wp) :: devk52 = 0. 117 REAL(wp) :: devk53 = 0.0794E-3 118 REAL(wp) :: devk54 = 0.09e-3 119 REAL(wp) :: devk55 = 0.054e-3 120 REAL(wp) :: devk56 = 0.3692E-3 121 REAL(wp) :: devk57 = 0.0427e-3 122 REAL(wp) :: devk58 = 0.09e-3 123 REAL(wp) :: devk59 = 0.0714e-3 124 REAL(wp) :: devk510 = 0.0 125 ! 126 ! General parameters 127 REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 128 REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp 129 130 ! Maximum number of iterations for each method 131 INTEGER, PARAMETER :: jp_maxniter_atgen = 20 132 133 ! Bookkeeping variables for each method 134 ! - SOLVE_AT_GENERAL 135 INTEGER :: niter_atgen = jp_maxniter_atgen 98 136 99 137 !!* Substitution … … 115 153 !!--------------------------------------------------------------------- 116 154 INTEGER :: ji, jj, jk 117 REAL(wp) :: ztkel, zt , zt2, zsal , zsal2 , zbuf1 , zbuf2155 REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2 118 156 REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 119 157 REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2 120 158 REAL(wp) :: zsqrt, ztr , zlogt , zcek1, zc1, zplat 121 REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1 159 REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1, za2 122 160 REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw 161 REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 123 162 REAL(wp) :: zst , zft , zcks , zckf , zaksp1 163 REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total 164 124 165 !!--------------------------------------------------------------------- 125 166 ! 126 167 IF( nn_timing == 1 ) CALL timing_start('p4z_che') 168 ! 169 ! Computation of chemical constants require practical salinity 170 ! Thus, when TEOS08 is used, absolute salinity is converted to 171 ! practical salinity 172 ! ------------------------------------------------------------- 173 IF (nn_eos == -1) THEN 174 salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504 175 ELSE 176 salinprac(:,:,:) = tsn(:,:,:,jp_sal) 177 ENDIF 178 127 179 ! 128 180 ! Computations of chemical constants require in situ temperature … … 135 187 DO ji = 1, jpi 136 188 zpres = fsdept(ji,jj,jk) / 1000. 137 za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * ( tsn(ji,jj,jk,jp_sal) - 35.0) )189 za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) ) 138 190 za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 139 191 tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 … … 151 203 ztkel = tempis(ji,jj,1) + 273.15 152 204 zt = ztkel * 0.01 153 zt2 = zt * zt 154 zsal = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 155 zsal2 = zsal * zsal 156 zlogt = LOG( zt ) 205 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 157 206 ! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 158 207 ! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 159 208 zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel & 160 209 & + 0.0047036e-4*ztkel**2) 161 ! ! SET SOLUBILITIES OF O2 AND CO2 162 chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 210 chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 163 211 chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 164 212 chemc(ji,jj,3) = 57.7 - 0.118*ztkel … … 176 224 DO ji = 1, jpi 177 225 ztkel = tempis(ji,jj,jk) + 273.15 178 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.226 zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35. 179 227 zsal2 = zsal * zsal 180 228 ztgg = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature … … 183 231 ztgg4 = ztgg3 * ztgg 184 232 ztgg5 = ztgg4 * ztgg 185 zoxy = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5 & 186 + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) + ox10 * zsal2 233 234 zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 & 235 & + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 & 236 & - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) & 237 & - 3.11680e-7 * zsal2 187 238 chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm) 188 239 END DO … … 209 260 ! SET ABSOLUTE TEMPERATURE 210 261 ztkel = tempis(ji,jj,jk) + 273.15 211 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.262 zsal = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 212 263 zsqrt = SQRT( zsal ) 213 264 zsal15 = zsqrt * zsal … … 220 271 221 272 ! CHLORINITY (WOOSTER ET AL., 1969) 222 zcl = zsal * salchl273 zcl = zsal / 1.80655 223 274 224 275 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 225 zst = st1 * zcl * st2276 zst = 0.14 * zcl /96.062 226 277 227 278 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 228 zft = ft1 * zcl * ft2279 zft = 0.000067 * zcl /18.9984 229 280 230 281 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) … … 234 285 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 235 286 & + LOG(1.0 - 0.001005 * zsal)) 236 !237 aphscale(ji,jj,jk) = ( 1. + zst / zcks )238 287 239 288 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) … … 249 298 & * zlogt + 0.053105*zsqrt*ztkel 250 299 251 252 300 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 253 301 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale … … 257 305 - 0.01781*zsal + 0.0001122*zsal*zsal) 258 306 259 ! PKW (H2O) (DICKSON AND RILEY, 1979) 260 zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt & 261 & + (118.67*ztr - 5.977 + 1.0495 * zlogt) & 262 & * zsqrt - 0.01615 * zsal 307 ! PKW (H2O) (MILLERO, 1995) from composite data 308 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 309 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 310 311 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 312 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 313 & + (-106.736*ztr + 0.69171) * zsqrt & 314 & + (-0.65643*ztr - 0.01844) * zsal 315 316 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 317 & + (-160.340*ztr + 1.3566)*zsqrt & 318 & + (0.37335*ztr - 0.05778)*zsal 319 320 zck3p = -3070.75*ztr - 18.126 & 321 & + (17.27039*ztr + 2.81197) * zsqrt & 322 & + (-44.99486*ztr - 0.09984) * zsal 323 324 ! CONSTANT FOR SILICATE, MILLERO (1995) 325 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 326 & + (-458.79*ztr + 3.5913) * zisqrt & 327 & + (188.74*ztr - 1.5998) * zis & 328 & + (-12.1652*ztr + 0.07871) * zis2 & 329 & + LOG(1.0 - 0.001005*zsal) 263 330 264 331 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER … … 268 335 & - 0.07711*zsal + 0.0041249*zsal15 269 336 337 ! CONVERT FROM DIFFERENT PH SCALES 338 total2free = 1.0/(1.0 + zst/zcks) 339 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 340 total2SWS = total2free * free2SWS 341 SWS2total = 1.0 / total2SWS 342 270 343 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 271 zak1 = 10**(zck1) 272 zak2 = 10**(zck2) 273 zakb = EXP( zckb )344 zak1 = 10**(zck1) * total2SWS 345 zak2 = 10**(zck2) * total2SWS 346 zakb = EXP( zckb ) * total2SWS 274 347 zakw = EXP( zckw ) 275 348 zaksp1 = 10**(zaksp0) 349 zak1p = exp( zck1p ) 350 zak2p = exp( zck2p ) 351 zak3p = exp( zck3p ) 352 zaksi = exp( zcksi ) 353 zckf = zckf * total2SWS 276 354 277 355 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) … … 285 363 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 286 364 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 287 zcpexp = zpres / (rgas*ztkel)288 zcpexp2 = zpres * z pres/(rgas*ztkel)365 zcpexp = zpres / (rgas*ztkel) 366 zcpexp2 = zpres * zcpexp 289 367 290 368 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE … … 292 370 ! (CF. BROECKER ET AL., 1982) 293 371 294 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 372 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 373 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 374 ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 375 376 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 295 377 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 296 ak 13(ji,jj,jk) = zak1* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )378 ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 297 379 298 380 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 299 381 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 300 ak 23(ji,jj,jk) = zak2* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )382 akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 301 383 302 384 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 303 385 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 304 ak b3(ji,jj,jk) = zakb* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )386 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 305 387 306 388 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 307 389 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 308 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 309 390 aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 391 392 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 393 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 394 akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 395 396 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 397 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 398 ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 399 400 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 401 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 402 ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 403 404 zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 405 zbuf2 = 0.5 * ( devk49 + devk59 * ztc ) 406 ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 407 408 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 409 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 410 aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 411 412 ! CONVERT FROM DIFFERENT PH SCALES 413 total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 414 free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 415 total2SWS = total2free * free2SWS 416 SWS2total = 1.0 / total2SWS 417 418 ! Convert to total scale 419 ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total 420 ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total 421 akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total 422 akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total 423 ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 424 ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 425 ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 426 aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 427 akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free 310 428 311 429 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 312 430 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 313 431 ! (P. 1285) AND BERNER (1976) 314 zbuf1 = - ( devk1 5 + devk25 * ztc + devk35* ztc * ztc )315 zbuf2 = 0.5 * ( devk4 5 + devk55* ztc )432 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 433 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 316 434 aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 317 435 318 ! TOTAL BORATE CONCENTR. [MOLES/L] 319 borat(ji,jj,jk) = bor1 * zcl * bor2 436 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 437 borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 438 sulfat(ji,jj,jk) = zst 439 fluorid(ji,jj,jk) = zft 320 440 321 441 ! Iron and SIO3 saturation concentration from ... 322 442 sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 323 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) ) 324 443 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 444 445 ! Liu and Millero (1999) only valid 5 - 50 degC 446 ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16 447 fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254.0/ztkel1)) 448 fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320.0/ztkel1) ) 449 fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(zis**0.5)) - (1996.0/ztkel1) ) 450 fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086.0/ztkel1) ) 451 fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980.0/ztkel1) ) 325 452 END DO 326 453 END DO … … 331 458 END SUBROUTINE p4z_che 332 459 460 SUBROUTINE p4z_che_ahini( p_hini ) 461 !!--------------------------------------------------------------------- 462 !! *** ROUTINE ahini_for_at *** 463 !! 464 !! Subroutine returns the root for the 2nd order approximation of the 465 !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 466 !! polynomial) around the local minimum, if it exists. 467 !! Returns * 1E-03_wp if p_alkcb <= 0 468 !! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 469 !! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 470 !! and the 2nd order approximation does not have 471 !! a solution 472 !!--------------------------------------------------------------------- 473 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini 474 INTEGER :: ji, jj, jk 475 REAL(wp) :: zca1, zba1 476 REAL(wp) :: zd, zsqrtd, zhmin 477 REAL(wp) :: za2, za1, za0 478 REAL(wp) :: p_dictot, p_bortot, p_alkcb 479 480 IF( nn_timing == 1 ) CALL timing_start('p4z_che_ahini') 481 ! 482 DO jk = 1, jpk 483 DO jj = 1, jpj 484 DO ji = 1, jpi 485 p_alkcb = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 486 p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) 487 p_bortot = borat(ji,jj,jk) 488 IF (p_alkcb <= 0.) THEN 489 p_hini(ji,jj,jk) = 1.e-3 490 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 491 p_hini(ji,jj,jk) = 1.e-10_wp 492 ELSE 493 zca1 = p_dictot/( p_alkcb + rtrn ) 494 zba1 = p_bortot/ (p_alkcb + rtrn ) 495 ! Coefficients of the cubic polynomial 496 za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 497 za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1) & 498 & + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 499 za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 500 ! Taylor expansion around the minimum 501 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 502 ! for the minimum close to the root 503 504 IF(zd > 0.) THEN ! If the discriminant is positive 505 zsqrtd = SQRT(zd) 506 IF(za2 < 0) THEN 507 zhmin = (-za2 + zsqrtd)/3. 508 ELSE 509 zhmin = -za1/(za2 + zsqrtd) 510 ENDIF 511 p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 512 ELSE 513 p_hini(ji,jj,jk) = 1.e-7 514 ENDIF 515 ! 516 ENDIF 517 END DO 518 END DO 519 END DO 520 ! 521 IF( nn_timing == 1 ) CALL timing_stop('p4z_che_ahini') 522 ! 523 END SUBROUTINE p4z_che_ahini 524 525 !=============================================================================== 526 SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup ) 527 528 ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 529 ! contributions to total alkalinity (the infimum and the supremum), i.e 530 ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 531 532 ! Argument variables 533 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf 534 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 535 536 p_alknw_inf(:,:,:) = -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) & 537 & - fluorid(:,:,:) 538 p_alknw_sup(:,:,:) = (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) ) & 539 & * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 540 541 END SUBROUTINE anw_infsup 542 543 544 SUBROUTINE p4z_che_solve_hi( p_hini, zhi ) 545 546 ! Universal pH solver that converges from any given initial value, 547 ! determines upper an lower bounds for the solution if required 548 549 ! Argument variables 550 !-------------------- 551 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini 552 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi 553 554 ! Local variables 555 !----------------- 556 INTEGER :: ji, jj, jk, jn 557 REAL(wp) :: zh_ini, zh, zh_prev, zh_lnfactor 558 REAL(wp) :: zdelta, zh_delta 559 REAL(wp) :: zeqn, zdeqndh, zalka 560 REAL(wp) :: aphscale 561 REAL(wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 562 REAL(wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 563 REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 564 REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 565 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 566 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 567 REAL(wp) :: zalk_wat, zdalk_wat 568 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 569 LOGICAL :: l_exitnow 570 REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 572 573 IF( nn_timing == 1 ) CALL timing_start('p4z_che_solve_hi') 574 ! Allocate temporary workspace 575 CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 576 CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 577 578 CALL anw_infsup( zalknw_inf, zalknw_sup ) 579 580 rmask(:,:,:) = tmask(:,:,:) 581 zhi(:,:,:) = 0. 582 583 ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 584 DO jk = 1, jpk 585 DO jj = 1, jpj 586 DO ji = 1, jpi 587 IF (rmask(ji,jj,jk) == 1.) THEN 588 p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 589 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 590 zh_ini = p_hini(ji,jj,jk) 591 592 zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 593 594 IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 595 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 596 ELSE 597 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 598 ENDIF 599 600 zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 601 602 IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 603 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 604 ELSE 605 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 606 ENDIF 607 608 zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 609 ENDIF 610 END DO 611 END DO 612 END DO 613 614 zeqn_absmin(:,:,:) = HUGE(1._wp) 615 616 DO jn = 1, jp_maxniter_atgen 617 DO jk = 1, jpk 618 DO jj = 1, jpj 619 DO ji = 1, jpi 620 IF (rmask(ji,jj,jk) == 1.) THEN 621 zfact = rhop(ji,jj,jk) / 1000. + rtrn 622 p_alktot = trb(ji,jj,jk,jptal) / zfact 623 zdic = trb(ji,jj,jk,jpdic) / zfact 624 zbot = borat(ji,jj,jk) 625 zpt = trb(ji,jj,jk,jppo4) / zfact * po4r 626 zsit = trb(ji,jj,jk,jpsil) / zfact 627 zst = sulfat (ji,jj,jk) 628 zft = fluorid(ji,jj,jk) 629 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 630 zh = zhi(ji,jj,jk) 631 zh_prev = zh 632 633 ! H2CO3 - HCO3 - CO3 : n=2, m=0 634 znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 635 zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 636 zalk_dic = zdic * (znumer_dic/zdenom_dic) 637 zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh & 638 *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 639 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 640 641 642 ! B(OH)3 - B(OH)4 : n=1, m=0 643 znumer_bor = akb3(ji,jj,jk) 644 zdenom_bor = akb3(ji,jj,jk) + zh 645 zalk_bor = zbot * (znumer_bor/zdenom_bor) 646 zdnumer_bor = akb3(ji,jj,jk) 647 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 648 649 650 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 651 znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 652 & + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 653 zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 654 & + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 655 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 656 zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 657 & + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 658 & + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 659 & + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) & 660 & + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 661 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 662 663 ! H4SiO4 - H3SiO4 : n=1, m=0 664 znumer_sil = aksi3(ji,jj,jk) 665 zdenom_sil = aksi3(ji,jj,jk) + zh 666 zalk_sil = zsit * (znumer_sil/zdenom_sil) 667 zdnumer_sil = aksi3(ji,jj,jk) 668 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 669 670 ! HSO4 - SO4 : n=1, m=1 671 aphscale = 1.0 + zst/aks3(ji,jj,jk) 672 znumer_so4 = aks3(ji,jj,jk) * aphscale 673 zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 674 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 675 zdnumer_so4 = aks3(ji,jj,jk) 676 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 677 678 ! HF - F : n=1, m=1 679 znumer_flu = akf3(ji,jj,jk) 680 zdenom_flu = akf3(ji,jj,jk) + zh 681 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 682 zdnumer_flu = akf3(ji,jj,jk) 683 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 684 685 ! H2O - OH 686 aphscale = 1.0 + zst/aks3(ji,jj,jk) 687 zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale 688 zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 689 690 ! CALCULATE [ALK]([CO3--], [HCO3-]) 691 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 692 & + zalk_so4 + zalk_flu & 693 & + zalk_wat - p_alktot 694 695 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 696 & + zalk_so4 + zalk_flu + zalk_wat) 697 698 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 699 & + zdalk_so4 + zdalk_flu + zdalk_wat 700 701 ! Adapt bracketing interval 702 IF(zeqn > 0._wp) THEN 703 zh_min(ji,jj,jk) = zh_prev 704 ELSEIF(zeqn < 0._wp) THEN 705 zh_max(ji,jj,jk) = zh_prev 706 ENDIF 707 708 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 709 ! if the function evaluation at the current point is 710 ! not decreasing faster than with a bisection step (at least linearly) 711 ! in absolute value take one bisection step on [ph_min, ph_max] 712 ! ph_new = (ph_min + ph_max)/2d0 713 ! 714 ! In terms of [H]_new: 715 ! [H]_new = 10**(-ph_new) 716 ! = 10**(-(ph_min + ph_max)/2d0) 717 ! = SQRT(10**(-(ph_min + phmax))) 718 ! = SQRT(zh_max * zh_min) 719 zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 720 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 721 ELSE 722 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 723 ! = -zdeqndh * LOG(10) * [H] 724 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 725 ! 726 ! pH_new = pH_old + \deltapH 727 ! 728 ! [H]_new = 10**(-pH_new) 729 ! = 10**(-pH_old - \Delta pH) 730 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 731 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 732 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 733 734 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 735 736 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 737 zh = zh_prev*EXP(zh_lnfactor) 738 ELSE 739 zh_delta = zh_lnfactor*zh_prev 740 zh = zh_prev + zh_delta 741 ENDIF 742 743 IF( zh < zh_min(ji,jj,jk) ) THEN 744 ! if [H]_new < [H]_min 745 ! i.e., if ph_new > ph_max then 746 ! take one bisection step on [ph_prev, ph_max] 747 ! ph_new = (ph_prev + ph_max)/2d0 748 ! In terms of [H]_new: 749 ! [H]_new = 10**(-ph_new) 750 ! = 10**(-(ph_prev + ph_max)/2d0) 751 ! = SQRT(10**(-(ph_prev + phmax))) 752 ! = SQRT([H]_old*10**(-ph_max)) 753 ! = SQRT([H]_old * zh_min) 754 zh = SQRT(zh_prev * zh_min(ji,jj,jk)) 755 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 756 ENDIF 757 758 IF( zh > zh_max(ji,jj,jk) ) THEN 759 ! if [H]_new > [H]_max 760 ! i.e., if ph_new < ph_min, then 761 ! take one bisection step on [ph_min, ph_prev] 762 ! ph_new = (ph_prev + ph_min)/2d0 763 ! In terms of [H]_new: 764 ! [H]_new = 10**(-ph_new) 765 ! = 10**(-(ph_prev + ph_min)/2d0) 766 ! = SQRT(10**(-(ph_prev + ph_min))) 767 ! = SQRT([H]_old*10**(-ph_min)) 768 ! = SQRT([H]_old * zhmax) 769 zh = SQRT(zh_prev * zh_max(ji,jj,jk)) 770 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 771 ENDIF 772 ENDIF 773 774 zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 775 776 ! Stop iterations once |\delta{[H]}/[H]| < rdel 777 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 778 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 779 780 ! Alternatively: 781 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 782 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 783 ! < 1/LOG(10) * rdel 784 785 ! Hence |zeqn/(zdeqndh*zh)| < rdel 786 787 ! rdel <-- pp_rdel_ah_target 788 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 789 790 IF(l_exitnow) THEN 791 rmask(ji,jj,jk) = 0. 792 ENDIF 793 794 zhi(ji,jj,jk) = zh 795 796 IF(jn >= jp_maxniter_atgen) THEN 797 zhi(ji,jj,jk) = -1._wp 798 ENDIF 799 800 ENDIF 801 END DO 802 END DO 803 END DO 804 END DO 805 ! 806 CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 807 CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 808 809 810 IF( nn_timing == 1 ) CALL timing_stop('p4z_che_solve_hi') 811 812 813 END SUBROUTINE p4z_che_solve_hi 333 814 334 815 INTEGER FUNCTION p4z_che_alloc() … … 336 817 !! *** ROUTINE p4z_che_alloc *** 337 818 !!---------------------------------------------------------------------- 338 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), & 339 & tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 819 INTEGER :: ierr(3) ! Local variables 820 !!---------------------------------------------------------------------- 821 822 ierr(:) = 0 823 824 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) ) 825 826 ALLOCATE( akb3(jpi,jpj,jpk) , tempis(jpi, jpj, jpk), & 827 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , & 828 & aks3(jpi,jpj,jpk) , akf3(jpi,jpj,jpk) , & 829 & ak1p3(jpi,jpj,jpk) , ak2p3(jpi,jpj,jpk) , & 830 & ak3p3(jpi,jpj,jpk) , aksi3(jpi,jpj,jpk) , & 831 & fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) , & 832 & salinprac(jpi,jpj,jpk), STAT=ierr(2) ) 833 834 ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) ) 835 836 !* Variable for chemistry of the CO2 cycle 837 p4z_che_alloc = MAXVAL( ierr ) 340 838 ! 341 839 IF( p4z_che_alloc /= 0 ) CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') … … 355 853 356 854 !!====================================================================== 357 END MODULE p4zche855 END MODULE p4zche -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
r6943 r7607 87 87 REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact 88 88 REAL(wp) :: zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 89 REAL(wp) :: zph, z ah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co289 REAL(wp) :: zph, zdic, zsch_o2, zsch_co2 90 90 REAL(wp) :: zyr_dec, zdco2dt 91 91 CHARACTER (len=25) :: charout … … 122 122 #endif 123 123 124 DO jm = 1, 10 125 !CDIR NOVERRCHK 126 DO jj = 1, jpj 127 !CDIR NOVERRCHK 128 DO ji = 1, jpi 129 130 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 131 zbot = borat(ji,jj,1) 132 zfact = rhop(ji,jj,1) / 1000. + rtrn 133 zdic = trb(ji,jj,1,jpdic) / zfact 134 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 135 zalka = trb(ji,jj,1,jptal) / zfact 136 137 ! CALCULATE [ALK]([CO3--], [HCO3-]) 138 zalk = zalka - ( akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1) & 139 & + zbot / ( 1.+ zph / akb3(ji,jj,1) ) ) 140 141 ! CALCULATE [H+] AND [H2CO3] 142 zah2 = SQRT( (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1) & 143 & / ak13(ji,jj,1) ) * ( 2.* zdic - zalk ) ) 144 zah2 = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 ) 145 zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact 146 hi(ji,jj,1) = zah2 * zfact 147 END DO 124 DO jj = 1, jpj 125 DO ji = 1, jpi 126 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 127 zfact = rhop(ji,jj,1) / 1000. + rtrn 128 zdic = trb(ji,jj,1,jpdic) 129 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 130 ! CALCULATE [H2CO3] 131 zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 148 132 END DO 149 133 END DO 150 151 134 152 135 ! -------------- … … 254 237 ENDIF 255 238 ! 239 #if defined key_cpl_carbon_cycle 240 ! change units for carbon cycle coupling 241 oce_co2(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r ! in molC/m2/s 242 #endif 243 ! 256 244 CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 257 245 ! -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
r6943 r7607 23 23 USE sms_pisces ! PISCES Source Minus Sink variables 24 24 USE prtctl_trc ! print control for debugging 25 USE p4zche ! Chemical model 25 26 USE iom ! I/O manager 26 27 … … 60 61 ! 61 62 INTEGER, INTENT(in) :: kt, knt ! ocean time step 62 INTEGER :: ji, jj, jk, jn 63 REAL(wp) :: zalk, zdic, zph, zah2 64 REAL(wp) :: zdispot, zfact, zcalcon, zalka, zaldi 63 INTEGER :: ji, jj, jk 64 REAL(wp) :: zdispot, zfact, zcalcon 65 65 REAL(wp) :: zomegaca, zexcess, zexcess0 66 66 CHARACTER (len=25) :: charout 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss, zhinit, zhi 68 68 !!--------------------------------------------------------------------- 69 69 ! 70 70 IF( nn_timing == 1 ) CALL timing_start('p4z_lys') 71 71 ! 72 CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )72 CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss, zhinit, zhi ) 73 73 ! 74 74 zco3 (:,:,:) = 0. 75 75 zcaldiss(:,:,:) = 0. 76 zhinit(:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn ) 76 77 ! ------------------------------------------- 77 78 ! COMPUTE [CO3--] and [H+] CONCENTRATIONS 78 79 ! ------------------------------------------- 79 80 DO jn = 1, 5 ! BEGIN OF ITERATION 81 ! 82 !CDIR NOVERRCHK 83 DO jk = 1, jpkm1 84 !CDIR NOVERRCHK 85 DO jj = 1, jpj 86 !CDIR NOVERRCHK 87 DO ji = 1, jpi 88 zfact = rhop(ji,jj,jk) / 1000. + rtrn 89 zph = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+] 90 zdic = trb(ji,jj,jk,jpdic) / zfact 91 zalka = trb(ji,jj,jk,jptal) / zfact 92 ! CALCULATE [ALK]([CO3--], [HCO3-]) 93 zalk = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn ) & 94 & + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 95 ! CALCULATE [H+] and [CO3--] 96 zaldi = zdic - zalk 97 zah2 = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) ) 98 zah2 = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 ) 99 ! 100 zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact 101 hi(ji,jj,jk) = zah2 * zfact 102 END DO 80 81 CALL p4z_che_solve_hi( zhinit, zhi ) 82 83 DO jk = 1, jpkm1 84 DO jj = 1, jpj 85 DO ji = 1, jpi 86 zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2 & 87 & + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn ) 88 hi(ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000. 103 89 END DO 104 90 END DO 105 ! 106 END DO 91 END DO 107 92 108 93 ! --------------------------------------------------------- … … 138 123 ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION 139 124 zcaldiss(ji,jj,jk) = zdispot * rfact2 / rmtss ! calcite dissolution 140 zco3(ji,jj,jk) = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)141 125 ! 142 126 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk) … … 167 151 ENDIF 168 152 ! 169 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zc o3sat, zcaldiss)153 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat ) 170 154 ! 171 155 IF( nn_timing == 1 ) CALL timing_stop('p4z_lys') … … 186 170 !! 187 171 !!---------------------------------------------------------------------- 188 INTEGER :: ji, jj, jk189 172 INTEGER :: ios ! Local integer output status for namelist read 190 REAL(wp) :: zcaralk, zbicarb, zco3191 REAL(wp) :: ztmas, ztmas1192 193 173 NAMELIST/nampiscal/ kdca, nca 194 174 !!---------------------------------------------------------------------- -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
r6420 r7607 85 85 CALL p4z_che ! initialize the chemical constants 86 86 ! 87 IF( .NOT. ln_rsttr ) THEN ; CALL p4z_ ph_ini ! set PH at kt=nit00087 IF( .NOT. ln_rsttr ) THEN ; CALL p4z_che_ahini( hi ) ! set PH at kt=nit000 88 88 ELSE ; CALL p4z_rst( nittrc000, 'READ' ) !* read or initialize all required fields 89 89 ENDIF … … 310 310 END SUBROUTINE p4z_sms_init 311 311 312 SUBROUTINE p4z_ph_ini313 !!---------------------------------------------------------------------314 !! *** ROUTINE p4z_ini_ph ***315 !!316 !! ** Purpose : Initialization of chemical variables of the carbon cycle317 !!---------------------------------------------------------------------318 INTEGER :: ji, jj, jk319 REAL(wp) :: zcaralk, zbicarb, zco3320 REAL(wp) :: ztmas, ztmas1321 !!---------------------------------------------------------------------322 323 ! Set PH from total alkalinity, borat (???), akb3 (???) and ak23 (???)324 ! --------------------------------------------------------325 DO jk = 1, jpk326 DO jj = 1, jpj327 DO ji = 1, jpi328 ztmas = tmask(ji,jj,jk)329 ztmas1 = 1. - tmask(ji,jj,jk)330 zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / ( 1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) ) )331 zco3 = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1332 zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )333 hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1334 END DO335 END DO336 END DO337 !338 END SUBROUTINE p4z_ph_ini339 340 312 SUBROUTINE p4z_rst( kt, cdrw ) 341 313 !!--------------------------------------------------------------------- … … 350 322 INTEGER , INTENT(in) :: kt ! ocean time-step 351 323 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 352 !353 INTEGER :: ji, jj, jk354 REAL(wp) :: zcaralk, zbicarb, zco3355 REAL(wp) :: ztmas, ztmas1356 324 !!--------------------------------------------------------------------- 357 325 … … 365 333 CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:) ) 366 334 ELSE 367 ! hi(:,:,:) = 1.e-9 368 CALL p4z_ph_ini 335 CALL p4z_che_ahini( hi ) 369 336 ENDIF 370 337 CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
r6287 r7607 93 93 94 94 !!* Variable for chemistry of the CO2 cycle 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ???96 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak13 !: ??? 97 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak23 !: ??? 98 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksp !: ??? 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ???100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ???101 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hi !: ??? 102 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: excess !: ??? … … 153 150 154 151 !* Variable for chemistry of the CO2 cycle 155 ALLOCATE( ak b3(jpi,jpj,jpk) , ak13 (jpi,jpj,jpk) ,&152 ALLOCATE( ak13 (jpi,jpj,jpk) , & 156 153 & ak23(jpi,jpj,jpk) , aksp (jpi,jpj,jpk) , & 157 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , &158 154 & hi (jpi,jpj,jpk) , excess(jpi,jpj,jpk) , & 159 155 & aphscale(jpi,jpj,jpk), STAT=ierr(4) ) -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
r7522 r7607 131 131 INTEGER, INTENT(in) :: kt 132 132 INTEGER :: jn 133 REAL(wp) :: zkt 133 REAL(wp) :: zkt, zrec 134 134 CHARACTER(len=1) :: cl1 ! 1 character 135 135 CHARACTER(len=2) :: cl2 ! 2 characters … … 153 153 ! 154 154 ! !* Restart: read in restart file 155 IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 156 iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. & 157 iom_varid( numrtr, 'ktdcy' , ldstop = .FALSE. ) > 0 ) THEN 158 CALL iom_get( numrtr, 'ktdcy', zkt ) ! A mean of qsr 155 IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 & 156 & .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 & 157 & .AND. iom_varid( numrtr, 'ktdcy' , ldstop = .FALSE. ) > 0 & 158 & .AND. iom_varid( numrtr, 'nrdcy' , ldstop = .FALSE. ) > 0 ) THEN 159 160 CALL iom_get( numrtr, 'ktdcy', zkt ) 159 161 rsecfst = INT( zkt ) * rdttrc(1) 160 162 IF(lwp) WRITE(numout,*) 'trc_qsr_mean: qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s ' 161 163 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean ) ! A mean of qsr 162 DO jn = 1, nb_rec_per_day 163 IF( jn <= 9 ) THEN 164 WRITE(cl1,'(i1)') jn 165 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) ) ! A mean of qsr 166 ELSE 167 WRITE(cl2,'(i2.2)') jn 168 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) ) ! A mean of qsr 169 ENDIF 170 ENDDO 164 CALL iom_get( numrtr, 'nrdcy', zrec ) ! Number of record per days 165 IF( INT( zrec ) == nb_rec_per_day ) THEN 166 DO jn = 1, nb_rec_per_day 167 IF( jn <= 9 ) THEN 168 WRITE(cl1,'(i1)') jn 169 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) ) ! A mean of qsr 170 ELSE 171 WRITE(cl2,'(i2.2)') jn 172 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) ) ! A mean of qsr 173 ENDIF 174 ENDDO 175 ELSE 176 DO jn = 1, nb_rec_per_day 177 qsr_arr(:,:,jn) = qsr_mean(:,:) 178 ENDDO 179 ENDIF 171 180 ELSE !* no restart: set from nit000 values 172 181 IF(lwp) WRITE(numout,*) 'trc_qsr_mean: qsr_mean set to nit000 values' … … 185 194 llnew = ( rseclast - rsecfst ) .ge. rdt_sampl ! new shortwave to store 186 195 IF( llnew ) THEN 187 IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, &196 IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 188 197 & ' time = ', rseclast/3600.,'hours ' 189 198 rsecfst = rseclast … … 199 208 IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file kt =', kt 200 209 IF(lwp) WRITE(numout,*) '~~~~~~~' 201 zkt = REAL( kt, wp ) 202 CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 210 zkt = REAL( kt, wp ) 211 zrec = REAL( nb_rec_per_day, wp ) 212 CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 213 CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec ) 203 214 DO jn = 1, nb_rec_per_day 204 215 IF( jn <= 9 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.