Changeset 7162 for branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
- Timestamp:
- 2016-11-01T14:23:51+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
r7068 r7162 5 5 !!====================================================================== 6 6 !! History : 3.5 ! 2012-07 (O. Aumont, A. Tagliabue, C. Ethe) Original code 7 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 7 8 !!---------------------------------------------------------------------- 8 9 !! p4z_fechem : Compute remineralization/scavenging of iron … … 13 14 USE trc ! passive tracers common variables 14 15 USE sms_pisces ! PISCES Source Minus Sink variables 15 USE p4zopt ! optical model16 16 USE p4zche ! chemical model 17 17 USE p4zsbc ! Boundary conditions from sediments … … 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 voelker 28 LOGICAL :: ln_ligvar !: boolean for variable ligand concentration following Tagliabue and voelker 29 REAL(wp), PUBLIC :: xlam1 !: scavenging rate of Iron 30 REAL(wp), PUBLIC :: xlamdust !: scavenging rate of Iron by dust 31 REAL(wp), PUBLIC :: ligand !: ligand concentration in the ocean 32 33 !!gm Not DOCTOR norm !!! 27 !! * Shared module variables 28 LOGICAL :: ln_fechem !: boolean for complex iron chemistry following Tagliabue and voelker 29 LOGICAL :: ln_ligvar !: boolean for variable ligand concentration following Tagliabue and voelker 30 LOGICAL :: ln_fecolloid !: boolean for variable colloidal fraction 31 REAL(wp), PUBLIC :: xlam1 !: scavenging rate of Iron 32 REAL(wp), PUBLIC :: xlamdust !: scavenging rate of Iron by dust 33 REAL(wp), PUBLIC :: ligand !: ligand concentration in the ocean 34 REAL(wp), PUBLIC :: kfep !: rate constant for nanoparticle formation 35 34 36 REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth 35 37 … … 54 56 !! and one particulate form (ln_fechem) 55 57 !!--------------------------------------------------------------------- 56 INTEGER, INTENT(in) :: kt, knt ! ocean time step57 !58 INTEGER :: ji, jj, jk, jic59 CHARACTER (len=25) :: charout58 ! 59 INTEGER, INTENT(in) :: kt, knt ! ocean time step 60 ! 61 INTEGER :: ji, jj, jk, jic, jn 60 62 REAL(wp) :: zdep, zlam1a, zlam1b, zlamfac 61 REAL(wp) :: zkeq, zfeequi, zfesatur, zfecoll 63 REAL(wp) :: zkeq, zfeequi, zfesatur, zfecoll, fe3sol 62 64 REAL(wp) :: zdenom1, zscave, zaggdfea, zaggdfeb, zcoag 63 65 REAL(wp) :: ztrc, zdust 64 REAL(wp) :: zdenom , zdenom265 REAL(wp) , POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig66 REAL(wp) , POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP66 REAL(wp) :: zdenom2 67 REAL(wp) :: zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2 68 REAL(wp) :: zrum, zcodel, zargu, zlight 67 69 REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand 68 70 REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2 69 71 REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2 70 REAL(wp) :: ztfe, zoxy 72 REAL(wp) :: ztfe, zoxy, zhplus 73 REAL(wp) :: zaggliga, zaggligb 74 REAL(wp) :: dissol, zligco 75 CHARACTER (len=25) :: charout 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP 78 REAL(wp), POINTER, DIMENSION(:,: ) :: zstrn, zstrn2 71 79 !!--------------------------------------------------------------------- 72 80 ! 73 81 IF( nn_timing == 1 ) CALL timing_start('p4z_fechem') 74 82 ! 75 CALL wrk_alloc( jpi,jpj,jpk, zFe3, zFeL1, zTL1, ztotlig ) 83 ! Allocate temporary workspace 84 CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip ) 76 85 zFe3 (:,:,:) = 0. 77 86 zFeL1(:,:,:) = 0. 78 87 zTL1 (:,:,:) = 0. 79 88 IF( ln_fechem ) THEN 80 CALL wrk_alloc( jpi,jpj,jpk, zFe2, zFeL2, zTL2, zFeP ) 89 CALL wrk_alloc( jpi, jpj, zstrn, zstrn2 ) 90 CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP ) 81 91 zFe2 (:,:,:) = 0. 82 92 zFeL2(:,:,:) = 0. … … 92 102 ztotlig(:,:,:) = MIN( ztotlig(:,:,:), 10. ) 93 103 ELSE 94 ztotlig(:,:,:) = ligand * 1E9 104 IF( ln_ligand ) THEN ; ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9 105 ELSE ; ztotlig(:,:,:) = ligand * 1E9 106 ENDIF 95 107 ENDIF 96 108 97 109 IF( ln_fechem ) THEN 110 ! compute the day length depending on latitude and the day 111 zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp ) 112 zcodel = ASIN( SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp ) ) 113 114 ! day length in hours 115 zstrn(:,:) = 0. 116 DO jj = 1, jpj 117 DO ji = 1, jpi 118 zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad ) 119 zargu = MAX( -1., MIN( 1., zargu ) ) 120 zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 121 END DO 122 END DO 123 124 ! Maximum light intensity 125 zstrn2(:,:) = zstrn(:,:) / 24. 126 WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 127 zstrn(:,:) = 24. / zstrn(:,:) 128 98 129 ! ------------------------------------------------------------ 99 130 ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009) … … 101 132 ! Chemistry is supposed to be fast enough to be at equilibrium 102 133 ! ------------------------------------------------------------ 103 DO jk = 1, jpkm1 134 DO jn = 1, 2 135 DO jk = 1, jpkm1 104 136 DO jj = 1, jpj 105 137 DO ji = 1, jpi 138 zlight = etot(ji,jj,jk) * zstrn(ji,jj) * float(2-jn) 139 zzstrn2 = zstrn2(ji,jj) * float(2-jn) + (1. - zstrn2(ji,jj) ) * float(jn-1) 106 140 ! Calculate ligand concentrations : assume 2/3rd of excess goes to 107 141 ! strong ligands (L1) and 1/3rd to weak ligands (L2) … … 110 144 zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand 111 145 ! ionic strength from Millero et al. 1987 112 zionic = 19.9201 * tsn(ji,jj,jk,jp_sal) / ( 1000. - 1.00488 * tsn(ji,jj,jk,jp_sal) + rtrn )113 146 zph = -LOG10( MAX( hi(ji,jj,jk), rtrn) ) 114 zoxy = trb(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )147 zoxy = trb(ji,jj,jk,jpoxy) 115 148 ! Fe2+ oxydation rate from Santana-Casiano et al. (2005) 116 zkox = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( t sn(ji,jj,jk,jp_tem) + 273.15 ) &117 & - 0.04406 * SQRT( tsn(ji,jj,jk,jp_sal) ) - 0.002847 * tsn(ji,jj,jk,jp_sal)149 zkox = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 ) & 150 & - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk) 118 151 zkox = ( 10.** zkox ) * spd 119 152 zkox = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn ) 120 153 ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006) 121 zkph2 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk) + 2. ))154 zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj)) 122 155 zkph1 = zkph2 / 5. 123 156 ! pass the dfe concentration from PISCES … … 165 198 ENDIF 166 199 ! solve for the other Fe species 167 z Fe3(ji,jj,jk) = MAX( 0., zxs )168 z Fep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk)/ kpr ) )200 zzFe3 = MAX( 0., zxs ) 201 zzFep = MAX( 0., ( ks * zzFe3 / kpr ) ) 169 202 zkappa2 = ( kb2 + zkph2 ) / kl2 170 zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) ) 171 zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) ) 172 zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) ) 203 zzFeL2 = MAX( 0., ( zzFe3 * zTL2(ji,jj,jk) ) / ( zkappa2 + zzFe3 ) ) 204 zzFeL1 = MAX( 0., ( ztfe / zb - za / zb * zzFe3 - zc / zb * zzFeL2 ) ) 205 zzFe2 = MAX( 0., ( ( zkph1 * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) ) 206 zFe3(ji,jj,jk) = zFe3(ji,jj,jk) + zzFe3 * zzstrn2 207 zFe2(ji,jj,jk) = zFe2(ji,jj,jk) + zzFe2 * zzstrn2 208 zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2 209 zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2 210 zFeP(ji,jj,jk) = zFeP(ji,jj,jk) + zzFeP * zzstrn2 173 211 END DO 174 212 END DO 213 END DO 175 214 END DO 176 215 ELSE … … 198 237 ! 199 238 ENDIF 200 ! 239 201 240 zdust = 0. ! if no dust available 202 !203 241 DO jk = 1, jpkm1 204 242 DO jj = 1, jpj … … 212 250 zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9 213 251 ELSE 214 zfeequi = zFe3(ji,jj,jk) * 1E-9 215 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 252 IF (ln_fecolloid) THEN 253 zfeequi = zFe3(ji,jj,jk) * 1E-9 254 zhplus = max( rtrn, hi(ji,jj,jk) ) 255 fe3sol = fesol(ji,jj,jk,1) * ( fesol(ji,jj,jk,2) * zhplus**2 & 256 & + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4) & 257 & + fesol(ji,jj,jk,5) / zhplus ) 258 zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) ) 259 ELSE 260 zfeequi = zFe3(ji,jj,jk) * 1E-9 261 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 262 fe3sol = 0. 263 kfep = 0. 264 ENDIF 216 265 ENDIF 217 218 ztrc = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 266 ztrc = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 219 267 IF( ln_dust ) zdust = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s 220 268 zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc … … 232 280 zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. ) 233 281 zlamfac = MIN( 1. , zlamfac ) 234 !!gm very small BUG : it is unlikely but possible that gdept_n = 0 .....235 282 zdep = MIN( 1., 1000. / gdept_n(ji,jj,jk) ) 236 283 zlam1b = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) ) … … 242 289 ! ---------------------------------------------------------------- 243 290 zlam1a = ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk) & 244 & + ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc))291 & + ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) 245 292 zaggdfea = zlam1a * xstep * zfecoll 293 ! 246 294 zlam1b = 3.53E3 * trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) 247 295 zaggdfeb = zlam1b * xstep * zfecoll 248 296 ! 249 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag 297 ! precipitation of Fe3+, creation of nanoparticles 298 precip(ji,jj,jk) = MAX( 0., ( zfeequi - fe3sol ) ) * kfep * xstep 299 ! 300 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb & 301 & - zcoag - precip(ji,jj,jk) 250 302 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea 251 303 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb 304 IF( ln_ligand ) THEN 305 zligco = MAX( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) ) 306 zaggliga = zlam1a * xstep * zligco 307 zaggligb = zlam1b * xstep * zligco 308 tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk) 309 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb 310 ENDIF 252 311 END DO 253 312 END DO … … 256 315 ! Define the bioavailable fraction of iron 257 316 ! ---------------------------------------- 258 IF( ln_fechem ) THEN 259 biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 ) 260 ELSE 261 biron(:,:,:) = trb(:,:,:,jpfer) 317 IF( ln_fechem ) THEN ; biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 ) 318 ELSE ; biron(:,:,:) = trb(:,:,:,jpfer) 319 ENDIF 320 321 IF( ln_ligand .AND. .NOT.ln_fechem) THEN 322 plig(:,:,:) = MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trn(:,:,:,jpfer) +rtrn ) ) ) 323 plig(:,:,:) = MAX( 0. , plig(:,:,:) ) 262 324 ENDIF 263 325 264 326 ! Output of some diagnostics variables 265 327 ! --------------------------------- 266 IF( lk_iomput .AND. knt == nrdttrc ) THEN 328 IF( lk_iomput ) THEN 329 IF( knt == nrdttrc ) THEN 267 330 IF( iom_use("Fe3") ) CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+ 268 331 IF( iom_use("FeL1") ) CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1 … … 276 339 IF( iom_use("TL2") ) CALL iom_put("TL2" , zTL2 (:,:,:) * tmask(:,:,:) ) ! TL2 277 340 ENDIF 341 ENDIF 278 342 ENDIF 279 343 … … 284 348 ENDIF 285 349 ! 286 CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig ) 287 IF( ln_fechem ) CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP ) 350 CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip ) 351 IF( ln_fechem ) THEN 352 CALL wrk_dealloc( jpi, jpj, zstrn, zstrn2 ) 353 CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP ) 354 ENDIF 288 355 ! 289 356 IF( nn_timing == 1 ) CALL timing_stop('p4z_fechem') … … 304 371 !! 305 372 !!---------------------------------------------------------------------- 306 NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand373 NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep 307 374 INTEGER :: ios ! Local integer output status for namelist read 308 375 … … 320 387 WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer' 321 388 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 322 WRITE(numout,*) ' enable complex iron chemistry scheme ln_fechem =', ln_fechem 323 WRITE(numout,*) ' variable concentration of ligand ln_ligvar =', ln_ligvar 324 WRITE(numout,*) ' scavenging rate of Iron xlam1 =', xlam1 325 WRITE(numout,*) ' scavenging rate of Iron by dust xlamdust =', xlamdust 326 WRITE(numout,*) ' ligand concentration in the ocean ligand =', ligand 389 WRITE(numout,*) ' enable complex iron chemistry scheme ln_fechem =', ln_fechem 390 WRITE(numout,*) ' variable concentration of ligand ln_ligvar =', ln_ligvar 391 WRITE(numout,*) ' Variable colloidal fraction of Fe3+ ln_fecolloid =', ln_fecolloid 392 WRITE(numout,*) ' scavenging rate of Iron xlam1 =', xlam1 393 WRITE(numout,*) ' scavenging rate of Iron by dust xlamdust =', xlamdust 394 WRITE(numout,*) ' ligand concentration in the ocean ligand =', ligand 395 WRITE(numout,*) ' rate constant for nanoparticle formation kfep =', kfep 327 396 ENDIF 328 397 ! … … 353 422 ! 354 423 END SUBROUTINE p4z_fechem_init 355 356 424 !!====================================================================== 357 425 END MODULE p4zfechem
Note: See TracChangeset
for help on using the changeset viewer.
