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 @ 7068

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

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