New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10368 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zfechem.F90 – NEMO

Ignore:
Timestamp:
2018-12-03T12:45:01+01:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zfechem.F90

    r10069 r10368  
    2727   LOGICAL          ::   ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker 
    2828   LOGICAL          ::   ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker 
    29    LOGICAL          ::   ln_fecolloid !: boolean for variable colloidal fraction 
    3029   REAL(wp), PUBLIC ::   xlam1        !: scavenging rate of Iron  
    3130   REAL(wp), PUBLIC ::   xlamdust     !: scavenging rate of Iron by dust  
     
    6867      REAL(wp) ::   za, zb, zc, zkappa1, zkappa2, za0, za1, za2 
    6968      REAL(wp) ::   zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2 
    70       REAL(wp) ::   ztfe, zoxy, zhplus 
     69      REAL(wp) ::   ztfe, zoxy, zhplus, zxlam 
    7170      REAL(wp) ::   zaggliga, zaggligb 
    7271      REAL(wp) ::   dissol, zligco 
     72      REAL(wp) :: zrfact2 
    7373      CHARACTER (len=25) :: charout 
    7474      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1 
     75      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d 
    7576      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zFeL2, zTL2, zFe2, zFeP 
    7677      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::   zstrn, zstrn2 
     
    129130         ! ------------------------------------------------------------ 
    130131         DO jn = 1, 2 
    131           DO jk = 1, jpkm1 
    132             DO jj = 1, jpj 
    133                DO ji = 1, jpi 
    134                   zlight = etot(ji,jj,jk) * zstrn(ji,jj) * REAL( 2-jn, wp ) 
    135                   zzstrn2 = zstrn2(ji,jj) * REAL( 2-jn, wp ) + (1. - zstrn2(ji,jj) ) * REAL( jn-1, wp ) 
    136                   ! Calculate ligand concentrations : assume 2/3rd of excess goes to 
    137                   ! strong ligands (L1) and 1/3rd to weak ligands (L2) 
    138                   ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9 
    139                   zTL1(ji,jj,jk) =                0.000001 + 0.67 * ztligand 
    140                   zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand 
    141                   ! ionic strength from Millero et al. 1987 
    142                   zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) ) 
    143                   zoxy   = trb(ji,jj,jk,jpoxy) 
    144                   ! Fe2+ oxydation rate from Santana-Casiano et al. (2005) 
    145                   zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 )  & 
    146                     &    - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk) 
    147                   zkox   = ( 10.** zkox ) * spd 
    148                   zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn ) 
    149                   ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006) 
    150                   zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj)) 
    151                   zkph1 = zkph2 / 5. 
    152                   ! pass the dfe concentration from PISCES 
    153                   ztfe = trb(ji,jj,jk,jpfer) * 1e9 
    154                   ! ---------------------------------------------------------- 
    155                   ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION 
    156                   ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom.  
    157                   ! ---------------------------------------------------------- 
    158                   ! calculate some parameters 
    159                   za = 1 + ks / kpr 
    160                   zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn ) 
    161                   zc = 1 + zkph2 / ( zkox + rtrn ) 
    162                   zkappa1 = ( kb1 + zkph1 + kth ) / kl1 
    163                   zkappa2 = ( kb2 + zkph2 ) / kl2 
    164                   za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za 
    165                   za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za & 
     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 & 
    166167                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za 
    167                   za0 = -zkappa1 * zkappa2 * ztfe / za 
    168                   zp  = za1 - za2 * za2 / 3. 
    169                   zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0 
    170                   zp3 = zp / 3. 
    171                   zq2 = zq / 2. 
    172                   zd  = zp3 * zp3 * zp3 + zq2 * zq2 
    173                   zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. ) 
    174                   ! compute the roots 
    175                   IF( zp > 0.) THEN 
    176                      ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) ) 
    177                      zphi =  zq / ( 2. * zr * zr * zr )  
    178                      zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log(x + sqrt(x^2+1)) 
    179                      zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3. 
    180                   ELSE 
    181                      IF( zd > 0. ) THEN 
    182                         zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) ) 
    183                         ! zphi = ACOSH( zfff ) 
    184                         zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log(x + sqrt(x^2-1)) 
    185                         zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3. 
     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. 
    186181                     ELSE 
    187                         zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) ) 
    188                         zphi = ACOS( zfff ) 
    189                         DO jic = 1, 3 
    190                            zfunc = -2 * zr * COS( zphi / 3. + 2. * REAL( jic - 1, wp ) * rpi / 3. ) - za2 / 3. 
    191                            IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc 
    192                         END DO 
     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 
    193195                     ENDIF 
    194                   ENDIF 
    195                   ! solve for the other Fe species 
    196                   zzFe3 = MAX( 0., zxs ) 
    197                   zzFep = MAX( 0., ( ks * zzFe3 / kpr ) ) 
    198                   zkappa2 = ( kb2 + zkph2 ) / kl2 
    199                   zzFeL2 = MAX( 0., ( zzFe3 * zTL2(ji,jj,jk) ) / ( zkappa2 + zzFe3 ) ) 
    200                   zzFeL1 = MAX( 0., ( ztfe / zb - za / zb * zzFe3 - zc / zb * zzFeL2 ) ) 
    201                   zzFe2  = MAX( 0., ( ( zkph1 * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) ) 
    202                   zFe3(ji,jj,jk)  = zFe3(ji,jj,jk)  + zzFe3 * zzstrn2 
    203                   zFe2(ji,jj,jk)  = zFe2(ji,jj,jk)  + zzFe2 * zzstrn2 
    204                   zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2 
    205                   zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2 
    206                   zFeP(ji,jj,jk)  = zFeP(ji,jj,jk)  + zzFeP * zzstrn2 
     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 
    207209               END DO 
    208210            END DO 
    209          END DO 
    210211         END DO 
    211212      ELSE 
     
    218219            DO jj = 1, jpj 
    219220               DO ji = 1, jpi 
    220                   zTL1(ji,jj,jk) = ztotlig(ji,jj,jk) 
    221                   zkeq           = fekeq(ji,jj,jk) 
    222                   zfesatur       = zTL1(ji,jj,jk) * 1E-9 
    223                   ztfe           = trb(ji,jj,jk,jpfer)  
     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)  
    224225                  ! Fe' is the root of a 2nd order polynom 
    225226                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               & 
    226                      &             + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       & 
    227                      &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq ) 
     227                     &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       & 
     228                     &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq ) 
    228229                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9 
    229230                  zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) ) 
     
    242243               ! Scavenging onto dust is also included as evidenced from the DUNE experiments. 
    243244               ! -------------------------------------------------------------------------------------- 
     245               zhplus  = max( rtrn, hi(ji,jj,jk) ) 
     246               fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  & 
     247               &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     & 
     248               &         + fesol(ji,jj,jk,5) / zhplus ) 
    244249               IF( ln_fechem ) THEN 
    245250                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9 
    246251                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9 
     252                  precip(ji,jj,jk) = 0.0 
    247253               ELSE 
    248254                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
    249                   IF (ln_fecolloid) THEN 
    250                      zhplus   = max( rtrn, hi(ji,jj,jk) ) 
    251                      fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  & 
    252                      &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     & 
    253                      &         + fesol(ji,jj,jk,5) / zhplus ) 
    254                      zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) ) 
    255                   ELSE 
    256                      zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 
    257                      fe3sol  = 0. 
    258                   ENDIF 
     255                  zhplus  = max( rtrn, hi(ji,jj,jk) ) 
     256                  fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  & 
     257                  &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     & 
     258                  &         + 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 
    259262               ENDIF 
    260263               ! 
    261264               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6  
    262                IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s 
    263                zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc 
     265               IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) & 
     266               &  * EXP( -gdept_n(ji,jj,jk) / 540. ) 
     267               IF (ln_ligand) THEN 
     268                  zxlam  = xlam1 * MAX( 1.E-3, EXP(-2 * etot(ji,jj,jk) / 10. ) * (1. - EXP(-2 * trb(ji,jj,jk,jpoxy) / 100.E-6 ) )) 
     269               ELSE 
     270                  zxlam  = xlam1 * 1.0 
     271               ENDIF 
     272               zlam1b = 3.e-5 + xlamdust * zdust + zxlam * ztrc 
    264273               zscave = zfeequi * zlam1b * xstep 
    265274 
     
    267276               ! to later allocate scavenged iron to the different organic pools 
    268277               ! --------------------------------------------------------- 
    269                zdenom1 = xlam1 * trb(ji,jj,jk,jppoc) / zlam1b 
    270                zdenom2 = xlam1 * trb(ji,jj,jk,jpgoc) / zlam1b 
     278               zdenom1 = zxlam * trb(ji,jj,jk,jppoc) / zlam1b 
     279               zdenom2 = zxlam * trb(ji,jj,jk,jpgoc) / zlam1b 
    271280 
    272281               !  Increased scavenging for very high iron concentrations found near the coasts  
     
    276285               zlamfac = MIN( 1.  , zlamfac ) 
    277286               zdep    = MIN( 1., 1000. / gdept_n(ji,jj,jk) ) 
    278                zlam1b  = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) ) 
    279                zcoag   = zfeequi * zlam1b * xstep + 1E-4 * ( 1. - zlamfac ) * zdep * xstep * trb(ji,jj,jk,jpfer) 
     287               zcoag   = 1E-4 * ( 1. - zlamfac ) * zdep * xstep * trb(ji,jj,jk,jpfer) 
    280288 
    281289               !  Compute the coagulation of colloidal iron. This parameterization  
     
    283291               !  It requires certainly some more work as it is very poorly constrained. 
    284292               !  ---------------------------------------------------------------- 
    285                zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    & 
    286                    &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) ) 
     293               zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    & 
     294                   &      + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) ) 
    287295               zaggdfea = zlam1a * xstep * zfecoll 
    288296               ! 
    289                zlam1b = 3.53E3 *  trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) 
     297               zlam1b   = 3.53E3 * trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) 
    290298               zaggdfeb = zlam1b * xstep * zfecoll 
    291                ! 
    292                ! precipitation of Fe3+, creation of nanoparticles 
    293                precip(ji,jj,jk) = MAX( 0., ( zfeequi - fe3sol ) ) * kfep * xstep 
    294299               ! 
    295300               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb & 
     
    297302               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea 
    298303               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb 
     304               zscav3d(ji,jj,jk)   = zscave 
     305               zcoll3d(ji,jj,jk)   = zaggdfea + zaggdfeb 
    299306               ! 
    300307            END DO 
     
    317324                  ! 
    318325                  zlam1b   = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) 
    319                   zligco   = MAX( ( 0.1 * trb(ji,jj,jk,jplgw) ), ( trb(ji,jj,jk,jplgw) - fe3sol ) ) 
     326                  zligco   = 0.5 * trn(ji,jj,jk,jplgw) 
    320327                  zaggliga = zlam1a * xstep * zligco 
    321328                  zaggligb = zlam1b * xstep * zligco 
    322329                  tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk) 
    323330                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb 
     331                  zlcoll3d(ji,jj,jk)  = zaggliga + zaggligb 
    324332               END DO 
    325333            END DO 
     
    328336         IF( .NOT.ln_fechem) THEN 
    329337            plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) ) 
    330             plig(:,:,:) =  MAX( 0. , plig(:,:,:) ) 
    331338         ENDIF 
    332339         ! 
     
    336343      IF( lk_iomput ) THEN 
    337344         IF( knt == nrdttrc ) THEN 
     345         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s 
    338346         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+ 
    339347         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1 
    340348         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1 
    341349         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL 
    342          IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron 
     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 ) 
    343354         IF( ln_fechem ) THEN 
    344355            IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+ 
     
    380391      INTEGER ::   ios   ! Local integer  
    381392      !! 
    382       NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep  
     393      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand, kfep  
    383394      !!---------------------------------------------------------------------- 
    384395      ! 
     
    401412         WRITE(numout,*) '      enable complex iron chemistry scheme      ln_fechem    =', ln_fechem 
    402413         WRITE(numout,*) '      variable concentration of ligand          ln_ligvar    =', ln_ligvar 
    403          WRITE(numout,*) '      Variable colloidal fraction of Fe3+       ln_fecolloid =', ln_fecolloid 
    404414         WRITE(numout,*) '      scavenging rate of Iron                   xlam1        =', xlam1 
    405415         WRITE(numout,*) '      scavenging rate of Iron by dust           xlamdust     =', xlamdust 
     
    407417         WRITE(numout,*) '      rate constant for nanoparticle formation  kfep         =', kfep 
    408418      ENDIF 
     419      !  
     420      IF (ln_ligand .AND. ln_fechem) CALL ctl_stop( 'STOP', 'p4z_fechem_init: ln_ligand and ln_fechem are incompatible') 
    409421      ! 
    410422      IF( ln_fechem ) THEN             ! set some constants used by the complexe chemistry scheme 
Note: See TracChangeset for help on using the changeset viewer.