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

source: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 @ 3461

Last change on this file since 3461 was 3461, checked in by cetlod, 12 years ago

branch:2012/dev_r3438_LOCEAN15_PISLOB, minor changes in p4zfechem.F90, see ticket #972
File size: 19.1 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   !!----------------------------------------------------------------------
8#if defined key_pisces
9   !!----------------------------------------------------------------------
10   !!   'key_top'       and                                      TOP models
11   !!   'key_pisces'                                       PISCES bio-model
12   !!----------------------------------------------------------------------
13   !!   p4z_fechem       :  Compute remineralization/scavenging of iron
14   !!   p4z_fechem_init  :  Initialisation of parameters for remineralisation
15   !!   p4z_fechem_alloc :  Allocate remineralisation variables
16   !!----------------------------------------------------------------------
17   USE oce_trc         !  shared variables between ocean and passive tracers
18   USE trc             !  passive tracers common variables
19   USE sms_pisces      !  PISCES Source Minus Sink variables
20   USE p4zopt          !  optical model
21   USE p4zche          !  chemical model
22   USE p4zsbc          !  Boundary conditions from sediments
23   USE prtctl_trc      !  print control for debugging
24   USE iom             !  I/O manager
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   p4z_fechem      ! called in p4zbio.F90
30   PUBLIC   p4z_fechem_init ! called in trcsms_pisces.F90
31
32   !! * Shared module variables
33   LOGICAL          ::  ln_fechem  = .FALSE.    !: boolean for complex iron chemistry following Tagliabue and voelker
34   LOGICAL          ::  ln_ligvar  = .FALSE.    !: boolean for variable ligand concentration following Tagliabue and voelker
35   REAL(wp), PUBLIC ::  xlam1      = 0.005_wp   !: scavenging rate of Iron
36   REAL(wp), PUBLIC ::  xlamdust   = 100.0_wp   !: scavenging rate of Iron by dust
37   REAL(wp), PUBLIC ::  ligand     = 0.6E-9_wp  !: ligand concentration in the ocean
38
39   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
40
41   !!* Substitution
42#  include "top_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
45   !! $Id: p4zrem.F90 3160 2011-11-20 14:27:18Z cetlod $
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE p4z_fechem( kt, jnt )
51      !!---------------------------------------------------------------------
52      !!                     ***  ROUTINE p4z_fechem  ***
53      !!
54      !! ** Purpose :   Compute remineralization/scavenging of iron
55      !!
56      !! ** Method  :   2 different chemistry models are available for iron
57      !!                (1) The simple chemistry model of Aumont and Bopp (2006)
58      !!                    based on one ligand and one inorganic form
59      !!                (2) The complex chemistry model of Tagliabue and
60      !!                    Voelker (2009) based on 2 ligands, 2 inorganic forms
61      !!                    and one particulate form (ln_fechem)
62      !!---------------------------------------------------------------------
63      !
64      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
65      !
66      INTEGER  ::   ji, jj, jk, jic
67      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
68      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll
69      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
70      REAL(wp) ::   ztrc, zdust
71#if ! defined key_kriest
72      REAL(wp) ::   zdenom, zdenom2
73#endif
74      REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig
75      REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
76      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
77      REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
78      REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
79      REAL(wp) :: ztfe, zoxy
80      REAL(wp) :: zstep
81      CHARACTER (len=25) :: charout
82      !!---------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('p4z_fechem')
85      !
86      ! Allocate temporary workspace
87                       CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
88      IF( ln_fechem )  CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
89
90      ! Total ligand concentration : Ligands can be chosen to be constant or variable
91      ! Parameterization from Tagliabue and Voelker (2011)
92      ! -------------------------------------------------
93      IF( ln_ligvar ) THEN
94         ztotlig(:,:,:) =  0.09 * trn(:,:,:,jpdoc) * 1E6 + ligand * 1E9
95         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
96      ELSE
97         ztotlig(:,:,:) = ligand * 1E9
98      ENDIF
99
100      IF( ln_fechem ) THEN
101         ! ------------------------------------------------------------
102         ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
103         ! This model is based on two ligands, Fe2+, Fe3+ and Fep
104         ! Chemistry is supposed to be fast enough to be at equilibrium
105         ! ------------------------------------------------------------
106         DO jk = 1, jpkm1
107            DO jj = 1, jpj
108               DO ji = 1, jpi
109                  ! Calculate ligand concentrations : assume 2/3rd of excess goes to
110                  ! strong ligands (L1) and 1/3rd to weak ligands (L2)
111                  ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
112                  zTL1(ji,jj,jk) =                0.000001 + 0.67 * ztligand
113                  zTL2(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
114                  ! ionic strength from Millero et al. 1987
115                  zionic = 19.9201 * tsn(ji,jj,jk,jp_sal) / ( 1000. - 1.00488 * tsn(ji,jj,jk,jp_sal) + rtrn )
116                  zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
117                  zoxy   = trn(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
118                  ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
119                  zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tsn(ji,jj,1,jp_tem) + 273.15 )  &
120                    &    - 0.04406 * SQRT( tsn(ji,jj,jk,jp_sal) ) - 0.002847 * tsn(ji,jj,jk,jp_sal)
121                  zkox   = ( 10.** zkox ) * spd
122                  zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
123                  ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
124                  zkph2 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk) + 2. ) )
125                  zkph1 = zkph2 / 5.
126                  ! pass the dfe concentration from PISCES
127                  ztfe = trn(ji,jj,jk,jpfer) * 1e9
128                  ! ----------------------------------------------------------
129                  ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
130                  ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom.
131                  ! ----------------------------------------------------------
132                  ! calculate some parameters
133                  za = 1 + ks / kpr
134                  zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn )
135                  zc = 1 + zkph2 / ( zkox + rtrn )
136                  zkappa1 = ( kb1 + zkph1 + kth ) / kl1
137                  zkappa2 = ( kb2 + zkph2 ) / kl2
138                  za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
139                  za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za &
140                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za
141                  za0 = -zkappa1 * zkappa2 * ztfe / za
142                  zp  = za1 - za2 * za2 / 3.
143                  zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0
144                  zp3 = zp / 3.
145                  zq2 = zq / 2.
146                  zd  = zp3 * zp3 * zp3 + zq2 * zq2
147                  zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. )
148                  ! compute the roots
149                  IF( zp > 0.) THEN
150                     ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) )
151                     zphi =  zq / ( 2. * zr * zr * zr ) 
152                     zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log(x + sqrt(x^2+1))
153                     zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3.
154                  ELSE
155                     IF( zd > 0. ) THEN
156                        zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) )
157                        ! zphi = ACOSH( zfff )
158                        zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log(x + sqrt(x^2-1))
159                        zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3.
160                     ELSE
161                        zfff = MIN( 1., zq / ( 2. * zr * zr * zr ) )
162                        zphi = ACOS( zfff )
163                        DO jic = 1, 3
164                           zfunc = -2 * zr * COS( zphi / 3. + 2. * FLOAT( jic - 1 ) * rpi / 3. ) - za2 / 3.
165                           IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc
166                        END DO
167                     ENDIF
168                  ENDIF
169                  ! solve for the other Fe species
170                  zFe3(ji,jj,jk) = MAX( 0., zxs ) 
171                  zFep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk) / kpr ) )
172                  zkappa2 = ( kb2 + zkph2 ) / kl2
173                  zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) )
174                  zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) )
175                  zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) )
176               END DO
177            END DO
178         END DO
179      ELSE
180         ! ------------------------------------------------------------
181         ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006)
182         ! This model is based on one ligand and Fe'
183         ! Chemistry is supposed to be fast enough to be at equilibrium
184         ! ------------------------------------------------------------
185         DO jk = 1, jpkm1
186            DO jj = 1, jpj
187               DO ji = 1, jpi
188                  zTL1(ji,jj,jk) = ztotlig(ji,jj,jk)
189                  zkeq           = fekeq(ji,jj,jk)
190                  zfesatur       = zTL1(ji,jj,jk) * 1E-9
191                  ztfe           = trn(ji,jj,jk,jpfer) 
192                  ! Fe' is the root of a 2nd order polynom
193                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
194                     &             + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
195                     &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
196                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
197                  zFeL1(ji,jj,jk) = MAX( 0., trn(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
198              END DO
199            END DO
200         END DO
201         !
202      ENDIF
203
204
205!CDIR NOVERRCHK
206      DO jk = 1, jpkm1
207!CDIR NOVERRCHK
208         DO jj = 1, jpj
209!CDIR NOVERRCHK
210            DO ji = 1, jpi
211               zstep = xstep
212# if defined key_degrad
213               zstep = zstep * facvol(ji,jj,jk)
214# endif
215               ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water.
216               ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
217               ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
218               ! --------------------------------------------------------------------------------------
219               IF( ln_fechem ) THEN
220                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) ) * 1E-9
221                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
222               ELSE
223                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
224                  zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
225               ENDIF
226#if defined key_kriest
227               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
228#else
229               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
230#endif
231               zdust  = dust(ji,jj) / ( wdust * 30.42 * 0.035 ) * tmask(ji,jj,jk)
232               zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc
233               zscave = zfeequi * zlam1b * zstep
234
235               ! Compute the different ratios for scavenging of iron
236               ! to later allocate scavenged iron to the different organic pools
237               ! ---------------------------------------------------------
238#if  defined key_kriest
239               zdenom1 = xlam1 * trn(ji,jj,jk,jppoc) / zlam1b 
240#else
241               zdenom1 = xlam1 * trn(ji,jj,jk,jppoc) / zlam1b
242               zdenom2 = xlam1 * trn(ji,jj,jk,jpgoc) / zlam1b
243#endif
244
245               !  Increased scavenging for very high iron concentrations found near the coasts
246               !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
247               !  -----------------------------------------------------------
248               zlam1b  = xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
249               zcoag   = zfeequi * zlam1b * zstep
250
251               !  Compute the coagulation of colloidal iron. This parameterization
252               !  could be thought as an equivalent of colloidal pumping.
253               !  It requires certainly some more work as it is very poorly constrained.
254               !  ----------------------------------------------------------------
255               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
256               zlamfac = MIN( 1.  , zlamfac )
257               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
258               zlam1a  = ( 0.369  * 0.3 * trn(ji,jj,jk,jpdoc) + 102.4  * trn(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
259                   &   + ( 114.   * 0.3 * trn(ji,jj,jk,jpdoc) + 5.09E3 * trn(ji,jj,jk,jppoc) )
260#if defined key_kriest
261               zlam1a   = zlam1a + 1E-4 * ( 1. - zlamfac ) * zdep
262               zaggdfea = zlam1a * zstep * zfecoll
263               zaggdfeb = 0.
264               !
265               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
266               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
267#else
268               zlam1b = 3.53E3 *   trn(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) + 1E-4 * ( 1. - zlamfac ) * zdep 
269               zaggdfea = zlam1a * zstep * zfecoll
270               zaggdfeb = zlam1b * zstep * zfecoll
271               !
272               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
273               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
274               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
275#endif
276            END DO
277         END DO
278      END DO
279      !
280      !  Define the bioavailable fraction of iron
281      !  ----------------------------------------
282      IF( ln_fechem ) THEN
283          biron(:,:,:) = MAX( 0., trn(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
284      ELSE
285          biron(:,:,:) = trn(:,:,:,jpfer) 
286      ENDIF
287
288      !  Output of some diagnostics variables
289      !     ---------------------------------
290      IF( ln_diatrc .AND. lk_iomput ) THEN
291         IF( jnt == nrdttrc ) THEN
292            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
293            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
294            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
295            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
296            CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
297            IF( ln_fechem ) THEN
298               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
299               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
300               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
301               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
302            ENDIF
303         ENDIF
304      ENDIF
305
306      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
307         WRITE(charout, FMT="('fechem')")
308         CALL prt_ctl_trc_info(charout)
309         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
310      ENDIF
311      !
312                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
313      IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
314      !
315      IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')
316      !
317   END SUBROUTINE p4z_fechem
318
319
320   SUBROUTINE p4z_fechem_init
321      !!----------------------------------------------------------------------
322      !!                  ***  ROUTINE p4z_fechem_init  ***
323      !!
324      !! ** Purpose :   Initialization of iron chemistry parameters
325      !!
326      !! ** Method  :   Read the nampisfer namelist and check the parameters
327      !!      called at the first timestep
328      !!
329      !! ** input   :   Namelist nampisfer
330      !!
331      !!----------------------------------------------------------------------
332      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand 
333
334      REWIND( numnatp )                     ! read numnatp
335      READ  ( numnatp, nampisfer )
336
337      IF(lwp) THEN                         ! control print
338         WRITE(numout,*) ' '
339         WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
340         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
341         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem =', ln_fechem
342         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar =', ln_ligvar
343         WRITE(numout,*) '    scavenging rate of Iron                   xlam1     =', xlam1
344         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust  =', xlamdust
345         WRITE(numout,*) '    ligand concentration in the ocean         ligand    =', ligand
346      ENDIF
347      !
348      IF( ln_fechem ) THEN
349         ! initialization of some constants used by the complexe chemistry scheme
350         ! ----------------------------------------------------------------------
351         spd = 3600. * 24.
352         con = 1.E9
353         ! LIGAND KINETICS (values from Witter et al. 2000)
354         ! Weak (L2) ligands
355         ! Phaeophytin
356         kl2 = 12.2E5  * spd / con
357         kb2 = 12.3E-6 * spd
358         ! Strong (L1) ligands
359         ! Saccharides
360         ! kl1 = 12.2E5  * spd / con
361         ! kb1 = 12.3E-6 * spd
362         ! DFOB-
363         kl1 = 19.6e5  * spd / con
364         kb1 = 1.5e-6  * spd
365         ! pcp and remin of Fe3p
366         ks  = 0.075
367         kpr = 0.05
368         ! thermal reduction of Fe3
369         kth = 0.0048 * 24.
370         !
371      ENDIF
372      !
373   END SUBROUTINE p4z_fechem_init
374
375#else
376   !!======================================================================
377   !!  Dummy module :                                   No PISCES bio-model
378   !!======================================================================
379CONTAINS
380   SUBROUTINE p4z_fechem                    ! Empty routine
381   END SUBROUTINE p4z_fechem
382#endif 
383
384   !!======================================================================
385END MODULE p4zfechem
Note: See TracBrowser for help on using the repository browser.