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.
p4zfechem.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 @ 7617

Last change on this file since 7617 was 7617, checked in by aumont, 7 years ago

update diagnostics + changes in quota code

File size: 24.1 KB
RevLine 
[3443]1MODULE p4zfechem
2   !!======================================================================
3   !!                         ***  MODULE p4zfechem  ***
4   !! TOP :   PISCES Compute iron chemistry and scavenging
5   !!======================================================================
[3461]6   !! History :   3.5  !  2012-07 (O. Aumont, A. Tagliabue, C. Ethe) Original code
[6453]7   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
[3443]8   !!----------------------------------------------------------------------
[6453]9#if defined key_pisces || defined key_pisces_quota
[3443]10   !!----------------------------------------------------------------------
11   !!   'key_top'       and                                      TOP models
[6453]12   !!   'key_pisces*'                                       PISCES bio-model
[3443]13   !!----------------------------------------------------------------------
14   !!   p4z_fechem       :  Compute remineralization/scavenging of iron
15   !!   p4z_fechem_init  :  Initialisation of parameters for remineralisation
16   !!   p4z_fechem_alloc :  Allocate remineralisation variables
17   !!----------------------------------------------------------------------
18   USE oce_trc         !  shared variables between ocean and passive tracers
19   USE trc             !  passive tracers common variables
20   USE sms_pisces      !  PISCES Source Minus Sink variables
21   USE p4zopt          !  optical model
22   USE p4zche          !  chemical model
23   USE p4zsbc          !  Boundary conditions from sediments
24   USE prtctl_trc      !  print control for debugging
25   USE iom             !  I/O manager
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   p4z_fechem      ! called in p4zbio.F90
31   PUBLIC   p4z_fechem_init ! called in trcsms_pisces.F90
32
33   !! * Shared module variables
[4147]34   LOGICAL          ::  ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
35   LOGICAL          ::  ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
[6453]36   LOGICAL          ::  ln_fecolloid !: boolean for variable colloidal fraction
[4147]37   REAL(wp), PUBLIC ::  xlam1        !: scavenging rate of Iron
38   REAL(wp), PUBLIC ::  xlamdust     !: scavenging rate of Iron by dust
39   REAL(wp), PUBLIC ::  ligand       !: ligand concentration in the ocean
[6453]40   REAL(wp), PUBLIC ::  kfep         !: rate constant for nanoparticle formation
[3443]41
[3446]42   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
[3443]43
44   !!* Substitution
45#  include "top_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
48   !! $Id: p4zrem.F90 3160 2011-11-20 14:27:18Z cetlod $
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
[6453]53   SUBROUTINE p4z_fechem( kt, jnt )
[3443]54      !!---------------------------------------------------------------------
55      !!                     ***  ROUTINE p4z_fechem  ***
56      !!
57      !! ** Purpose :   Compute remineralization/scavenging of iron
58      !!
59      !! ** Method  :   2 different chemistry models are available for iron
60      !!                (1) The simple chemistry model of Aumont and Bopp (2006)
61      !!                    based on one ligand and one inorganic form
62      !!                (2) The complex chemistry model of Tagliabue and
63      !!                    Voelker (2009) based on 2 ligands, 2 inorganic forms
64      !!                    and one particulate form (ln_fechem)
65      !!---------------------------------------------------------------------
66      !
[6453]67      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
[3443]68      !
[6453]69      INTEGER  ::   ji, jj, jk, jic, jn
[3446]70      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
[7617]71      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, zfecollc, fe3sol, fe3sol1
[3446]72      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
[3443]73      REAL(wp) ::   ztrc, zdust
74#if ! defined key_kriest
[6841]75      REAL(wp) ::   zdenom2
[3443]76#endif
[6453]77      REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip
[3443]78      REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
[6453]79      REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zstrn2
[7617]80      REAL(wp), POINTER, DIMENSION(:,:,:) :: zcoll3d, zscav3d
[6453]81      REAL(wp) ::   zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2
[6841]82      REAL(wp) ::   zrum, zcodel, zargu, zlight
[3446]83      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
[3443]84      REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
85      REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
[6453]86      REAL(wp) :: ztfe, zoxy, zhplus
[7617]87      REAL(wp) :: zstep, zrfact2
[6453]88#if defined key_ligand
89      REAL(wp) :: zaggliga, zaggligb
90      REAL(wp) :: dissol, zligco
[7617]91      REAL(wp), POINTER, DIMENSION(:,:,:) :: zlcoll3d
[6453]92#endif
[3443]93      CHARACTER (len=25) :: charout
94      !!---------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )  CALL timing_start('p4z_fechem')
97      !
98      ! Allocate temporary workspace
[6453]99      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
[7617]100      CALL wrk_alloc( jpi, jpj, jpk, zcoll3d, zscav3d )
101#if defined key_ligand
102      CALL wrk_alloc( jpi, jpj, jpk, zlcoll3d )
103#endif
[3531]104      zFe3 (:,:,:) = 0.
105      zFeL1(:,:,:) = 0.
106      zTL1 (:,:,:) = 0.
107      IF( ln_fechem ) THEN
[6453]108         CALL wrk_alloc( jpi, jpj,      zstrn, zstrn2 )
[3531]109         CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
110         zFe2 (:,:,:) = 0.
111         zFeL2(:,:,:) = 0.
112         zTL2 (:,:,:) = 0.
113         zFeP (:,:,:) = 0.
114      ENDIF
[3461]115
[3443]116      ! Total ligand concentration : Ligands can be chosen to be constant or variable
117      ! Parameterization from Tagliabue and Voelker (2011)
118      ! -------------------------------------------------
119      IF( ln_ligvar ) THEN
[5385]120         ztotlig(:,:,:) =  0.09 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9
[3446]121         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
[3443]122      ELSE
[6453]123#if defined key_ligand
124         ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9
125#else
[3443]126         ztotlig(:,:,:) = ligand * 1E9
[6453]127#endif
[3443]128      ENDIF
129
130      IF( ln_fechem ) THEN
[6453]131
132         ! compute the day length depending on latitude and the day
133         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
134         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
135
136         ! day length in hours
137         zstrn(:,:) = 0.
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
141               zargu = MAX( -1., MIN(  1., zargu ) )
142               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
143            END DO
144         END DO
145
146         ! Maximum light intensity
147         zstrn2(:,:) = zstrn(:,:) / 24.
148         WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
149         zstrn(:,:) = 24. / zstrn(:,:)
150
[3443]151         ! ------------------------------------------------------------
152         ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
153         ! This model is based on two ligands, Fe2+, Fe3+ and Fep
154         ! Chemistry is supposed to be fast enough to be at equilibrium
155         ! ------------------------------------------------------------
[6453]156         DO jn = 1, 2
[3531]157!CDIR NOVERRCHK
[3443]158         DO jk = 1, jpkm1
[3531]159!CDIR NOVERRCHK
[3443]160            DO jj = 1, jpj
[3531]161!CDIR NOVERRCHK
[3443]162               DO ji = 1, jpi
[6453]163                  zlight = etot(ji,jj,jk) * zstrn(ji,jj) * float(2-jn)
164                  zzstrn2 = zstrn2(ji,jj) * float(2-jn) + (1. - zstrn2(ji,jj) ) * float(jn-1)
[3461]165                  ! Calculate ligand concentrations : assume 2/3rd of excess goes to
166                  ! strong ligands (L1) and 1/3rd to weak ligands (L2)
[3446]167                  ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
[3461]168                  zTL1(ji,jj,jk) =                0.000001 + 0.67 * ztligand
169                  zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
[3443]170                  ! ionic strength from Millero et al. 1987
171                  zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
[6966]172                  zoxy   = trb(ji,jj,jk,jpoxy)
[3461]173                  ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
[6966]174                  zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 )  &
175                    &    - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk)
[3461]176                  zkox   = ( 10.** zkox ) * spd
177                  zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
[3443]178                  ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
[6841]179                  zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj))
[3461]180                  zkph1 = zkph2 / 5.
[3443]181                  ! pass the dfe concentration from PISCES
[5385]182                  ztfe = trb(ji,jj,jk,jpfer) * 1e9
[3443]183                  ! ----------------------------------------------------------
184                  ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
185                  ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom.
186                  ! ----------------------------------------------------------
187                  ! calculate some parameters
188                  za = 1 + ks / kpr
[3446]189                  zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn )
[3443]190                  zc = 1 + zkph2 / ( zkox + rtrn )
[3446]191                  zkappa1 = ( kb1 + zkph1 + kth ) / kl1
192                  zkappa2 = ( kb2 + zkph2 ) / kl2
[3443]193                  za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
194                  za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za &
195                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za
196                  za0 = -zkappa1 * zkappa2 * ztfe / za
197                  zp  = za1 - za2 * za2 / 3.
198                  zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0
199                  zp3 = zp / 3.
200                  zq2 = zq / 2.
201                  zd  = zp3 * zp3 * zp3 + zq2 * zq2
202                  zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. )
203                  ! compute the roots
204                  IF( zp > 0.) THEN
[3450]205                     ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) )
206                     zphi =  zq / ( 2. * zr * zr * zr ) 
[3461]207                     zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log(x + sqrt(x^2+1))
[3443]208                     zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3.
209                  ELSE
210                     IF( zd > 0. ) THEN
211                        zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) )
[3461]212                        ! zphi = ACOSH( zfff )
213                        zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log(x + sqrt(x^2-1))
[3443]214                        zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3.
215                     ELSE
[3461]216                        zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) )
[3456]217                        zphi = ACOS( zfff )
[3443]218                        DO jic = 1, 3
219                           zfunc = -2 * zr * COS( zphi / 3. + 2. * FLOAT( jic - 1 ) * rpi / 3. ) - za2 / 3.
220                           IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc
221                        END DO
222                     ENDIF
223                  ENDIF
224                  ! solve for the other Fe species
[6453]225                  zzFe3 = MAX( 0., zxs )
226                  zzFep = MAX( 0., ( ks * zzFe3 / kpr ) )
[3448]227                  zkappa2 = ( kb2 + zkph2 ) / kl2
[6453]228                  zzFeL2 = MAX( 0., ( zzFe3 * zTL2(ji,jj,jk) ) / ( zkappa2 + zzFe3 ) )
229                  zzFeL1 = MAX( 0., ( ztfe / zb - za / zb * zzFe3 - zc / zb * zzFeL2 ) )
230                  zzFe2  = MAX( 0., ( ( zkph1 * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) )
231                  zFe3(ji,jj,jk)  = zFe3(ji,jj,jk)  + zzFe3 * zzstrn2
232                  zFe2(ji,jj,jk)  = zFe2(ji,jj,jk)  + zzFe2 * zzstrn2
233                  zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2
234                  zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2
235                  zFeP(ji,jj,jk)  = zFeP(ji,jj,jk)  + zzFeP * zzstrn2
[3443]236               END DO
237            END DO
238         END DO
[6453]239         END DO
[3443]240      ELSE
241         ! ------------------------------------------------------------
242         ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006)
243         ! This model is based on one ligand and Fe'
244         ! Chemistry is supposed to be fast enough to be at equilibrium
245         ! ------------------------------------------------------------
[3531]246!CDIR NOVERRCHK
[3443]247         DO jk = 1, jpkm1
[3531]248!CDIR NOVERRCHK
[3443]249            DO jj = 1, jpj
[3531]250!CDIR NOVERRCHK
[3443]251               DO ji = 1, jpi
252                  zTL1(ji,jj,jk) = ztotlig(ji,jj,jk)
253                  zkeq           = fekeq(ji,jj,jk)
254                  zfesatur       = zTL1(ji,jj,jk) * 1E-9
[5385]255                  ztfe           = trb(ji,jj,jk,jpfer) 
[3443]256                  ! Fe' is the root of a 2nd order polynom
[3448]257                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
[3443]258                     &             + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
259                     &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
[3448]260                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
[5385]261                  zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
[3443]262              END DO
263            END DO
264         END DO
265         !
266      ENDIF
267
[3531]268      zdust = 0.         ! if no dust available
[3443]269!CDIR NOVERRCHK
270      DO jk = 1, jpkm1
271!CDIR NOVERRCHK
272         DO jj = 1, jpj
273!CDIR NOVERRCHK
274            DO ji = 1, jpi
275               zstep = xstep
276# if defined key_degrad
277               zstep = zstep * facvol(ji,jj,jk)
278# endif
279               ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water.
280               ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
281               ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
282               ! --------------------------------------------------------------------------------------
[3446]283               IF( ln_fechem ) THEN
[4521]284                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9
[3446]285                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
286               ELSE
[7180]287                  zfeequi = zFe3(ji,jj,jk) * 1E-9
[6453]288                  IF (ln_fecolloid) THEN
289                     zhplus   = max( rtrn, hi(ji,jj,jk) )
[7180]290                     fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
[6453]291                     &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
292                     &         + fesol(ji,jj,jk,5) / zhplus )
293                     zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) )
[7180]294#if defined key_ligand
295                     zligco  = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) )
296#endif
[6453]297                  ELSE
298                     zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
[7180]299#if defined key_ligand
300                     zligco  = 0.5 * trn(ji,jj,jk,jplgw)
301#endif
[6453]302                     fe3sol  = 0.
303                  ENDIF
[3446]304               ENDIF
[3443]305#if defined key_kriest
[5385]306               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
[3443]307#else
[5385]308               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
[3443]309#endif
[4800]310               IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s
[3443]311               zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc
312               zscave = zfeequi * zlam1b * zstep
313
314               ! Compute the different ratios for scavenging of iron
315               ! to later allocate scavenged iron to the different organic pools
316               ! ---------------------------------------------------------
[5385]317               zdenom1 = xlam1 * trb(ji,jj,jk,jppoc) / zlam1b
[3780]318#if ! defined key_kriest
[5385]319               zdenom2 = xlam1 * trb(ji,jj,jk,jpgoc) / zlam1b
[3443]320#endif
321
322               !  Increased scavenging for very high iron concentrations found near the coasts
323               !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
324               !  -----------------------------------------------------------
[3475]325               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
326               zlamfac = MIN( 1.  , zlamfac )
327               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
[5385]328               zlam1b  = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
329               zcoag   = zfeequi * zlam1b * zstep + 1E-4 * ( 1. - zlamfac ) * zdep * zstep * trb(ji,jj,jk,jpfer)
[3443]330
331               !  Compute the coagulation of colloidal iron. This parameterization
332               !  could be thought as an equivalent of colloidal pumping.
333               !  It requires certainly some more work as it is very poorly constrained.
334               !  ----------------------------------------------------------------
[7617]335               zfecollc = trb(ji,jj,jk,jpdoc) / ( 40.E-6 + trb(ji,jj,jk,jpdoc) ) * zfecoll
[5385]336               zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
[7180]337                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
[7617]338               zaggdfea = zlam1a * zstep * zfecollc
[6453]339               !
[3443]340#if defined key_kriest
[3446]341               zaggdfeb = 0.
[3443]342#else
[5385]343               zlam1b = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
[7617]344               zaggdfeb = zlam1b * zstep * zfecollc
[7180]345#endif
[6453]346               ! precipitation of Fe3+, creation of nanoparticles
347               precip(ji,jj,jk) = max( 0., (zfeequi - fe3sol) ) * kfep * zstep
348               !
349               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb &
350               &                     - zcoag - precip(ji,jj,jk)
[3446]351               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
[7180]352#if ! defined key_kriest
[3446]353               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
[3443]354#endif
[7617]355               zscav3d(ji,jj,jk)  = zscave
356               zcoll3d(ji,jj,jk)  = zaggdfea + zaggdfeb
[6453]357#if defined key_ligand
[7180]358               zaggliga = zlam1a * zstep * zligco
359#   if defined key_kriest
360               zaggligb = 0.
361#   else
362               zaggligb = zlam1b * zstep * zligco
363#   endif
[6453]364               tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk)
365               tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
[7617]366               zlcoll3d(ji,jj,jk) = zaggliga + zaggligb
[6453]367#endif
[3443]368            END DO
369         END DO
370      END DO
371      !
[3446]372      !  Define the bioavailable fraction of iron
373      !  ----------------------------------------
374      IF( ln_fechem ) THEN
[5385]375          biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
[3446]376      ELSE
[5385]377          biron(:,:,:) = trb(:,:,:,jpfer) 
[3446]378      ENDIF
[3443]379
[6453]380#if defined key_ligand
381      IF (.NOT. ln_fechem) THEN
382          plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trn(:,:,:,jpfer) +rtrn ) ) )
383          plig(:,:,:) = max(0. , plig(:,:,:) )
384      ENDIF
385#endif
386
[3443]387      !  Output of some diagnostics variables
388      !     ---------------------------------
[6453]389      IF( ln_diatrc .AND. lk_iomput ) THEN
390         IF( jnt == nrdttrc ) THEN
[7617]391            zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
[6453]392            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
393            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
394            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
395            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
[7617]396            CALL iom_put("Biron"  , biron  (:,:,:)  * 1e9 * tmask(:,:,:) )   ! biron
397            CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
398            CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
399#if defined key_ligand
400            CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
401#endif
[6453]402            IF( ln_fechem ) THEN
403               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
404               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
405               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
406               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
407            ENDIF
[3443]408         ENDIF
409      ENDIF
410
411      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
[3449]412         WRITE(charout, FMT="('fechem')")
[3443]413         CALL prt_ctl_trc_info(charout)
414         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
415      ENDIF
416      !
[7180]417      CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
[7617]418      CALL wrk_dealloc( jpi, jpj, jpk, zscav3d, zcoll3d )
419#if defined key_ligand
420      CALL wrk_dealloc( jpi, jpj, jpk, zlcoll3d )
421#endif
422
[7180]423      IF( ln_fechem )  THEN
424         CALL wrk_dealloc( jpi, jpj,      zstrn, zstrn2 )
425         CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
426      ENDIF
[3443]427      !
428      IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')
429      !
430   END SUBROUTINE p4z_fechem
431
432
433   SUBROUTINE p4z_fechem_init
434      !!----------------------------------------------------------------------
435      !!                  ***  ROUTINE p4z_fechem_init  ***
436      !!
437      !! ** Purpose :   Initialization of iron chemistry parameters
438      !!
439      !! ** Method  :   Read the nampisfer namelist and check the parameters
440      !!      called at the first timestep
441      !!
442      !! ** input   :   Namelist nampisfer
443      !!
444      !!----------------------------------------------------------------------
[6453]445      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep 
[4147]446      INTEGER :: ios                 ! Local integer output status for namelist read
[3443]447
[4147]448      REWIND( numnatp_ref )              ! Namelist nampisfer in reference namelist : Pisces iron chemistry
449      READ  ( numnatp_ref, nampisfer, IOSTAT = ios, ERR = 901)
450901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisfer in reference namelist', lwp )
[3443]451
[4147]452      REWIND( numnatp_cfg )              ! Namelist nampisfer in configuration namelist : Pisces iron chemistry
453      READ  ( numnatp_cfg, nampisfer, IOSTAT = ios, ERR = 902 )
454902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisfer in configuration namelist', lwp )
[4624]455      IF(lwm) WRITE ( numonp, nampisfer )
[4147]456
[3443]457      IF(lwp) THEN                         ! control print
458         WRITE(numout,*) ' '
459         WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
460         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
[6453]461         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem    =', ln_fechem
462         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar    =', ln_ligvar
463         WRITE(numout,*) '    Variable colloidal fraction of Fe3+       ln_fecolloid =', ln_fecolloid
464         WRITE(numout,*) '    scavenging rate of Iron                   xlam1        =', xlam1
465         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust     =', xlamdust
466         WRITE(numout,*) '    ligand concentration in the ocean         ligand       =', ligand
467         WRITE(numout,*) '    rate constant for nanoparticle formation  kfep         =', kfep
[3443]468      ENDIF
469      !
470      IF( ln_fechem ) THEN
471         ! initialization of some constants used by the complexe chemistry scheme
472         ! ----------------------------------------------------------------------
[3446]473         spd = 3600. * 24.
[3443]474         con = 1.E9
475         ! LIGAND KINETICS (values from Witter et al. 2000)
[3461]476         ! Weak (L2) ligands
477         ! Phaeophytin
478         kl2 = 12.2E5  * spd / con
479         kb2 = 12.3E-6 * spd
480         ! Strong (L1) ligands
481         ! Saccharides
482         ! kl1 = 12.2E5  * spd / con
483         ! kb1 = 12.3E-6 * spd
484         ! DFOB-
485         kl1 = 19.6e5  * spd / con
486         kb1 = 1.5e-6  * spd
[3443]487         ! pcp and remin of Fe3p
488         ks  = 0.075
489         kpr = 0.05
[3446]490         ! thermal reduction of Fe3
491         kth = 0.0048 * 24.
[3461]492         !
[3443]493      ENDIF
494      !
495   END SUBROUTINE p4z_fechem_init
496
497#else
498   !!======================================================================
499   !!  Dummy module :                                   No PISCES bio-model
500   !!======================================================================
501CONTAINS
502   SUBROUTINE p4z_fechem                    ! Empty routine
503   END SUBROUTINE p4z_fechem
504#endif 
505
506   !!======================================================================
507END MODULE p4zfechem
Note: See TracBrowser for help on using the repository browser.