Changeset 10401 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zfechem.F90
- Timestamp:
- 2018-12-17T16:57:34+01:00 (5 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/TOP/PISCES/P4Z/p4zfechem.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/TOP/PISCES/P4Z/p4zfechem.F90
r10362 r10401 25 25 PUBLIC p4z_fechem_init ! called in trcsms_pisces.F90 26 26 27 LOGICAL :: ln_fechem !: boolean for complex iron chemistry following Tagliabue and voelker28 27 LOGICAL :: ln_ligvar !: boolean for variable ligand concentration following Tagliabue and voelker 29 28 REAL(wp), PUBLIC :: xlam1 !: scavenging rate of Iron … … 32 31 REAL(wp), PUBLIC :: kfep !: rate constant for nanoparticle formation 33 32 34 REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth !!gm <<<== DOCTOR names SVP !!!35 36 33 !!---------------------------------------------------------------------- 37 34 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 47 44 !! ** Purpose : Compute remineralization/scavenging of iron 48 45 !! 49 !! ** Method : 2 different chemistry models are available for iron 50 !! (1) The simple chemistry model of Aumont and Bopp (2006) 51 !! based on one ligand and one inorganic form 52 !! (2) The complex chemistry model of Tagliabue and 53 !! Voelker (2009) based on 2 ligands, 2 inorganic forms 54 !! and one particulate form (ln_fechem) 46 !! ** Method : A simple chemistry model of iron from Aumont and Bopp (2006) 47 !! based on one ligand and one inorganic form 55 48 !!--------------------------------------------------------------------- 56 49 INTEGER, INTENT(in) :: kt, knt ! ocean time step … … 74 67 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTL1, zFe3, ztotlig, precip, zFeL1 75 68 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcoll3d, zscav3d, zlcoll3d 76 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zFeL2, zTL2, zFe2, zFeP77 REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zstrn, zstrn278 69 !!--------------------------------------------------------------------- 79 70 ! … … 83 74 zFeL1(:,:,:) = 0. 84 75 zTL1 (:,:,:) = 0. 85 IF( ln_fechem ) THEN86 ALLOCATE( zstrn(jpi,jpj), zstrn2(jpi,jpj) )87 ALLOCATE( zFe2(jpi,jpj,jpk), zFeL2(jpi,jpj,jpk), zTL2(jpi,jpj,jpk), zFeP(jpi,jpj,jpk) )88 zFe2 (:,:,:) = 0.89 zFeL2(:,:,:) = 0.90 zTL2 (:,:,:) = 0.91 zFeP (:,:,:) = 0.92 ENDIF93 76 94 77 ! Total ligand concentration : Ligands can be chosen to be constant or variable … … 104 87 ENDIF 105 88 106 IF( ln_fechem ) THEN 107 ! compute the day length depending on latitude and the day 108 zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp ) 109 zcodel = ASIN( SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp ) ) 110 111 ! day length in hours 112 zstrn(:,:) = 0. 89 ! ------------------------------------------------------------ 90 ! from Aumont and Bopp (2006) 91 ! This model is based on one ligand and Fe' 92 ! Chemistry is supposed to be fast enough to be at equilibrium 93 ! ------------------------------------------------------------ 94 DO jk = 1, jpkm1 113 95 DO jj = 1, jpj 114 96 DO ji = 1, jpi 115 zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad ) 116 zargu = MAX( -1., MIN( 1., zargu ) ) 117 zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 118 END DO 97 zTL1(ji,jj,jk) = ztotlig(ji,jj,jk) 98 zkeq = fekeq(ji,jj,jk) 99 zfesatur = zTL1(ji,jj,jk) * 1E-9 100 ztfe = trb(ji,jj,jk,jpfer) 101 ! Fe' is the root of a 2nd order polynom 102 zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe ) & 103 & + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2 & 104 & + 4. * ztfe * zkeq) ) / ( 2. * zkeq ) 105 zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9 106 zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) ) 107 END DO 119 108 END DO 120 121 ! Maximum light intensity 122 zstrn2(:,:) = zstrn(:,:) / 24. 123 WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 124 zstrn(:,:) = 24. / zstrn(:,:) 125 126 ! ------------------------------------------------------------ 127 ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009) 128 ! This model is based on two ligands, Fe2+, Fe3+ and Fep 129 ! Chemistry is supposed to be fast enough to be at equilibrium 130 ! ------------------------------------------------------------ 131 DO jn = 1, 2 132 DO jk = 1, jpkm1 133 DO jj = 1, jpj 134 DO ji = 1, jpi 135 zlight = etot(ji,jj,jk) * zstrn(ji,jj) * REAL( 2-jn, wp ) 136 zzstrn2 = zstrn2(ji,jj) * REAL( 2-jn, wp ) + (1. - zstrn2(ji,jj) ) * REAL( jn-1, wp ) 137 ! Calculate ligand concentrations : assume 2/3rd of excess goes to 138 ! strong ligands (L1) and 1/3rd to weak ligands (L2) 139 ztligand = ztotlig(ji,jj,jk) - ligand * 1E9 140 zTL1(ji,jj,jk) = 0.000001 + 0.67 * ztligand 141 zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand 142 ! ionic strength from Millero et al. 1987 143 zph = -LOG10( MAX( hi(ji,jj,jk), rtrn) ) 144 zoxy = trb(ji,jj,jk,jpoxy) 145 ! Fe2+ oxydation rate from Santana-Casiano et al. (2005) 146 zkox = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 ) & 147 & - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk) 148 zkox = ( 10.** zkox ) * spd 149 zkox = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn ) 150 ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006) 151 zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj)) 152 zkph1 = zkph2 / 5. 153 ! pass the dfe concentration from PISCES 154 ztfe = trb(ji,jj,jk,jpfer) * 1e9 155 ! ---------------------------------------------------------- 156 ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION 157 ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom. 158 ! ---------------------------------------------------------- 159 ! calculate some parameters 160 za = 1.0 + ks / kpr 161 zb = 1.0 + zkph2 / ( zkox ) 162 zc = 1.0 + ( zkph1 + kth ) / ( zkox ) 163 zkappa1 = ( kb1 + zkph1 + kth ) / kl1 164 zkappa2 = ( kb2 + zkph2 ) / kl2 165 za2 = zTL2(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za 166 za1 = zkappa1 * zTL2(ji,jj,jk) * zb / za + zkappa2 * zTL1(ji,jj,jk) * zc / za & 167 & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za 168 za0 = -zkappa1 * zkappa2 * ztfe / za 169 zp = za1 - za2 * za2 / 3. 170 zq = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0 171 zp3 = zp / 3. 172 zq2 = zq / 2. 173 zd = zp3 * zp3 * zp3 + zq2 * zq2 174 zr = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. ) 175 ! compute the roots 176 IF( zp > 0.) THEN 177 ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) ) 178 zphi = zq / ( 2. * zr * zr * zr ) 179 zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) ) ! asinh(x) = log(x + sqrt(x^2+1)) 180 zxs = -2. * zr * SINH( zphi / 3. ) - za1 / 3. 181 ELSE 182 IF( zd > 0. ) THEN 183 zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) ) 184 ! zphi = ACOSH( zfff ) 185 zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) ) ! acosh(x) = log(x + sqrt(x^2-1)) 186 zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3. 187 ELSE 188 zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) ) 189 zphi = ACOS( zfff ) 190 DO jic = 1, 3 191 zfunc = -2 * zr * COS( zphi / 3. + 2. * REAL( jic - 1, wp ) * rpi / 3. ) - za2 / 3. 192 IF( zfunc > 0. .AND. zfunc <= ztfe) zxs = zfunc 193 END DO 194 ENDIF 195 ENDIF 196 ! solve for the other Fe species 197 zzFe3 = MAX( 0., zxs ) 198 zzFep = MAX( 0., ( ks * zzFe3 / kpr ) ) 199 zzFeL1 = MAX( 0., ( zzFe3 * zTL1(ji,jj,jk) ) / ( zkappa1 + zzFe3 ) ) 200 zzFeL2 = (ztfe - za * zzFe3 - zc * zzFeL1 ) / zb 201 zzFe2 = MAX( 0., ( ( ( zkph1 + kth ) * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) ) 202 zzFep = ztfe - zzFe3 - zzFe2 - zzFeL1 - zzFeL2 203 zFe3(ji,jj,jk) = zFe3(ji,jj,jk) + zzFe3 * zzstrn2 204 zFe2(ji,jj,jk) = zFe2(ji,jj,jk) + zzFe2 * zzstrn2 205 zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2 206 zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2 207 zFeP(ji,jj,jk) = zFeP(ji,jj,jk) + zzFeP * zzstrn2 208 END DO 209 END DO 210 END DO 211 END DO 212 ELSE 213 ! ------------------------------------------------------------ 214 ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006) 215 ! This model is based on one ligand and Fe' 216 ! Chemistry is supposed to be fast enough to be at equilibrium 217 ! ------------------------------------------------------------ 218 DO jk = 1, jpkm1 219 DO jj = 1, jpj 220 DO ji = 1, jpi 221 zTL1(ji,jj,jk) = ztotlig(ji,jj,jk) 222 zkeq = fekeq(ji,jj,jk) 223 zfesatur = zTL1(ji,jj,jk) * 1E-9 224 ztfe = trb(ji,jj,jk,jpfer) 225 ! Fe' is the root of a 2nd order polynom 226 zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe ) & 227 & + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2 & 228 & + 4. * ztfe * zkeq) ) / ( 2. * zkeq ) 229 zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9 230 zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) ) 231 END DO 232 END DO 233 END DO 234 ! 235 ENDIF 109 END DO 110 ! 236 111 237 112 zdust = 0. ! if no dust available … … 247 122 & + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4) & 248 123 & + fesol(ji,jj,jk,5) / zhplus ) 249 IF( ln_fechem ) THEN 250 zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9 251 zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9 252 precip(ji,jj,jk) = 0.0 253 ELSE 254 zfeequi = zFe3(ji,jj,jk) * 1E-9 255 zhplus = max( rtrn, hi(ji,jj,jk) ) 256 fe3sol = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2 & 124 ! 125 zfeequi = zFe3(ji,jj,jk) * 1E-9 126 zhplus = max( rtrn, hi(ji,jj,jk) ) 127 fe3sol = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2 & 257 128 & + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4) & 258 129 & + fesol(ji,jj,jk,5) / zhplus ) 259 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 260 ! precipitation of Fe3+, creation of nanoparticles 261 precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep 262 ENDIF 130 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 131 ! precipitation of Fe3+, creation of nanoparticles 132 precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep 263 133 ! 264 134 ztrc = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 … … 311 181 ! Define the bioavailable fraction of iron 312 182 ! ---------------------------------------- 313 IF( ln_fechem ) THEN ; biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 ) 314 ELSE ; biron(:,:,:) = trb(:,:,:,jpfer) 315 ENDIF 183 biron(:,:,:) = trb(:,:,:,jpfer) 316 184 ! 317 185 IF( ln_ligand ) THEN … … 334 202 END DO 335 203 ! 336 IF( .NOT.ln_fechem) THEN 337 plig(:,:,:) = MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) ) 338 ENDIF 204 plig(:,:,:) = MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) ) 339 205 ! 340 206 ENDIF … … 343 209 IF( lk_iomput ) THEN 344 210 IF( knt == nrdttrc ) THEN 345 zrfact2 = 1.e3 * rfact2r ! conversion from mol/L/timestep into mol/m3/s 346 IF( iom_use("Fe3") ) CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+ 347 IF( iom_use("FeL1") ) CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1 348 IF( iom_use("TL1") ) CALL iom_put("TL1" , zTL1 (:,:,:) * tmask(:,:,:) ) ! TL1 349 IF( iom_use("Totlig") ) CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) ) ! TL 350 IF( iom_use("Biron") ) CALL iom_put("Biron" , biron (:,:,:) * 1e9 * tmask(:,:,:) ) ! biron 351 IF( iom_use("FESCAV") ) CALL iom_put("FESCAV" , zscav3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 352 IF( iom_use("FECOLL") ) CALL iom_put("FECOLL" , zcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 353 IF( iom_use("LGWCOLL")) CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 354 IF( ln_fechem ) THEN 355 IF( iom_use("Fe2") ) CALL iom_put("Fe2" , zFe2 (:,:,:) * tmask(:,:,:) ) ! Fe2+ 356 IF( iom_use("FeL2") ) CALL iom_put("FeL2" , zFeL2 (:,:,:) * tmask(:,:,:) ) ! FeL2 357 IF( iom_use("FeP") ) CALL iom_put("FeP" , zFeP (:,:,:) * tmask(:,:,:) ) ! FeP 358 IF( iom_use("TL2") ) CALL iom_put("TL2" , zTL2 (:,:,:) * tmask(:,:,:) ) ! TL2 359 ENDIF 211 zrfact2 = 1.e3 * rfact2r ! conversion from mol/L/timestep into mol/m3/s 212 IF( iom_use("Fe3") ) CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+ 213 IF( iom_use("FeL1") ) CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1 214 IF( iom_use("TL1") ) CALL iom_put("TL1" , zTL1 (:,:,:) * tmask(:,:,:) ) ! TL1 215 IF( iom_use("Totlig") ) CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) ) ! TL 216 IF( iom_use("Biron") ) CALL iom_put("Biron" , biron (:,:,:) * 1e9 * tmask(:,:,:) ) ! biron 217 IF( iom_use("FESCAV") ) CALL iom_put("FESCAV" , zscav3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 218 IF( iom_use("FECOLL") ) CALL iom_put("FECOLL" , zcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 219 IF( iom_use("LGWCOLL")) CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 ) 360 220 ENDIF 361 221 ENDIF … … 367 227 ENDIF 368 228 ! 369 IF( ln_fechem ) THEN370 DEALLOCATE( zstrn, zstrn2 )371 DEALLOCATE( zFe2, zFeL2, zTL2, zFeP )372 ENDIF373 !374 229 IF( ln_timing ) CALL timing_stop('p4z_fechem') 375 230 ! … … 391 246 INTEGER :: ios ! Local integer 392 247 !! 393 NAMELIST/nampisfer/ ln_ fechem, ln_ligvar, xlam1, xlamdust, ligand, kfep248 NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep 394 249 !!---------------------------------------------------------------------- 395 250 ! … … 410 265 IF(lwp) THEN ! control print 411 266 WRITE(numout,*) ' Namelist : nampisfer' 412 WRITE(numout,*) ' enable complex iron chemistry scheme ln_fechem =', ln_fechem413 267 WRITE(numout,*) ' variable concentration of ligand ln_ligvar =', ln_ligvar 414 268 WRITE(numout,*) ' scavenging rate of Iron xlam1 =', xlam1 … … 418 272 ENDIF 419 273 ! 420 IF (ln_ligand .AND. ln_fechem) CALL ctl_stop( 'STOP', 'p4z_fechem_init: ln_ligand and ln_fechem are incompatible')421 !422 IF( ln_fechem ) THEN ! set some constants used by the complexe chemistry scheme423 !424 spd = 3600. * 24.425 con = 1.E9426 ! LIGAND KINETICS (values from Witter et al. 2000)427 ! Weak (L2) ligands428 ! Phaeophytin429 kl2 = 12.2E5 * spd / con430 kb2 = 12.3E-6 * spd431 ! Strong (L1) ligands432 ! Saccharides433 ! kl1 = 12.2E5 * spd / con434 ! kb1 = 12.3E-6 * spd435 ! DFOB-436 kl1 = 19.6e5 * spd / con437 kb1 = 1.5e-6 * spd438 ! pcp and remin of Fe3p439 ks = 0.075440 kpr = 0.05441 ! thermal reduction of Fe3442 kth = 0.0048 * 24.443 !444 ENDIF445 !446 274 END SUBROUTINE p4z_fechem_init 447 275
Note: See TracChangeset
for help on using the changeset viewer.
