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

source: branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 @ 7041

Last change on this file since 7041 was 7041, checked in by cetlod, 8 years ago

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