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_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 @ 9461

Last change on this file since 9461 was 9461, checked in by aumont, 6 years ago

bug fixes in iron cycle

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