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