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

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

ranch:2012/dev_r3438_LOCEAN15_PISLOB : change ACOSH/ASINH intrinsic function with their mathematical formulae, see ticket #972

File size: 19.4 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, 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      ! Total ligand concentration : Ligands can be chosen to be constant or variable
90      ! Parameterization from Tagliabue and Voelker (2011)
91      ! -------------------------------------------------
92      IF( ln_ligvar ) THEN
93         ztotlig(:,:,:) =  0.09 * trn(:,:,:,jpdoc) * 1E6 + ligand * 1E9
94         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
95      ELSE
96         ztotlig(:,:,:) = ligand * 1E9
97      ENDIF
98
99      IF( ln_fechem ) THEN
100         ! ------------------------------------------------------------
101         ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
102         ! This model is based on two ligands, Fe2+, Fe3+ and Fep
103         ! Chemistry is supposed to be fast enough to be at equilibrium
104         ! ------------------------------------------------------------
105         DO jk = 1, jpkm1
106            DO jj = 1, jpj
107               DO ji = 1, jpi
108                  ! Calculate ligand concentrations : assume 2/3rd of excess goes to L2 and 1/3rd to L1
109                  ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
110                  zTL2(ji,jj,jk) =                0.000001 + 0.67 * ztligand
111                  zTL1(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
112                  ! ionic strength from Millero et al. 1987
113                  zionic = 19.9201 * tsn(ji,jj,jk,jp_sal) / ( 1000. - 1.00488 * tsn(ji,jj,jk,jp_sal) + rtrn )
114                  zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
115                  zoxy   = trn(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
116! from Trapp and Millero (2007) corrected
117!                 zkox=(10.**( &
118!                   &     ((-.02*(zionic**.5)))-(3841/(tn &
119!                   &   (ji,jj,jk)+273.15))-(3.76*pH)+(0.33* &
120!                   &   pH**2)+(1.15*log10(carb))+27.87))* &
121!                   &   max(oxy,1.e-6) &
122!                   kox = kox * (60*24)
123                ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
124                  zkox = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tsn(ji,jj,1,jp_tem) + 273.15 )  &
125                    &   - 0.04406 * SQRT( tsn(ji,jj,jk,jp_sal) ) - 0.002847 * tsn(ji,jj,jk,jp_sal)
126                  zkox = ( 10.** zkox ) * spd
127                  zkox = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
128                  ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
129                  ! kph1=etot(ji,jj,jk)*0.0000001672*2.857142857143*spd modified to enhance
130                  zkph1 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk) + 2. ) )
131                  zkph2 = zkph1 / 5.
132                  ! pass the dfe concentration from PISCES
133                  ztfe = trn(ji,jj,jk,jpfer) * 1e9
134                  ! ----------------------------------------------------------
135                  ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
136                  ! As shown in Tagliabue and Voelker (2009), Fe3+ is the root of a 3rd order polynom.
137                  ! ----------------------------------------------------------
138                  ! calculate some parameters
139                  za = 1 + ks / kpr
140                  zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn )
141                  zc = 1 + zkph2 / ( zkox + rtrn )
142                  zkappa1 = ( kb1 + zkph1 + kth ) / kl1
143                  zkappa2 = ( kb2 + zkph2 ) / kl2
144                  za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
145                  za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za &
146                      & + zkappa1 * zkappa2 - ( zkappa1 + zkappa2 ) * ztfe / za
147                  za0 = -zkappa1 * zkappa2 * ztfe / za
148                  zp  = za1 - za2 * za2 / 3.
149                  zq  = za2 * za2 * za2 * 2. / 27. - za2 * za1 / 3. + za0
150                  zp3 = zp / 3.
151                  zq2 = zq / 2.
152                  zd  = zp3 * zp3 * zp3 + zq2 * zq2
153                  zr  = zq / ABS( zq ) * SQRT( ABS( zp ) / 3. )
154                  ! compute the roots
155                  IF( zp > 0.) THEN
156                     ! zphi = ASINH( zq / ( 2. * zr * zr * zr ) )
157                     zphi =  zq / ( 2. * zr * zr * zr ) 
158                     zphi = LOG( zphi + SQRT( zphi * zphi + 1 ) )  ! asinh(x) = log( x + sqrt( x*x + 1 ) )
159                     zxs  = -2. * zr * SINH( zphi / 3. ) - za1 / 3.
160                  ELSE
161                     IF( zd > 0. ) THEN
162                        zfff = MAX( 1., zq / ( 2. * zr * zr * zr ) )
163                      !  zphi = ACOSH( zfff )
164                        zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log( x + sqrt( x*x - 1 ) )
165                        zxs = -2. * zr * COSH( zphi / 3. ) - za1 / 3.
166                     ELSE
167                        zfff = zq / ( 2. * zr * zr * zr )
168                        zfff = MIN( 1., zfff )
169                        !zphi = ACOSH( zfff )
170                        zphi = LOG( zfff + SQRT( zfff * zfff - 1 ) )  ! acosh(x) = log( x + sqrt( x*x - 1 ) )
171                        DO jic = 1, 3
172                           zfunc = -2 * zr * COS( zphi / 3. + 2. * FLOAT( jic - 1 ) * rpi / 3. ) - za2 / 3.
173                           IF( zfunc > 0. .AND. zfunc <= ztfe)  zxs = zfunc
174                        END DO
175                     ENDIF
176                  ENDIF
177                  ! solve for the other Fe species
178                  zFe3(ji,jj,jk) = MAX( 0., zxs ) 
179                  zFep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk) / kpr ) )
180                  zkappa2 = ( kb2 + zkph2 ) / kl2
181                  zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) )
182                  zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) )
183                  zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) )
184               END DO
185            END DO
186         END DO
187      ELSE
188         ! ------------------------------------------------------------
189         ! OLD FE CHEMISTRY ROUTINE from Aumont and Bopp (2006)
190         ! This model is based on one ligand and Fe'
191         ! Chemistry is supposed to be fast enough to be at equilibrium
192         ! ------------------------------------------------------------
193         DO jk = 1, jpkm1
194            DO jj = 1, jpj
195               DO ji = 1, jpi
196                  zTL1(ji,jj,jk) = ztotlig(ji,jj,jk)
197                  zkeq           = fekeq(ji,jj,jk)
198                  zfesatur       = zTL1(ji,jj,jk) * 1E-9
199                  ztfe           = trn(ji,jj,jk,jpfer) 
200                  ! Fe' is the root of a 2nd order polynom
201                  zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
202                     &             + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
203                     &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
204                  zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
205                  zFeL1(ji,jj,jk) = MAX( 0., trn(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
206              END DO
207            END DO
208         END DO
209         !
210      ENDIF
211
212
213!CDIR NOVERRCHK
214      DO jk = 1, jpkm1
215!CDIR NOVERRCHK
216         DO jj = 1, jpj
217!CDIR NOVERRCHK
218            DO ji = 1, jpi
219               zstep = xstep
220# if defined key_degrad
221               zstep = zstep * facvol(ji,jj,jk)
222# endif
223               ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water.
224               ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
225               ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
226               ! --------------------------------------------------------------------------------------
227               IF( ln_fechem ) THEN
228                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) ) * 1E-9
229                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
230               ELSE
231                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
232                  zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
233               ENDIF
234#if defined key_kriest
235               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
236#else
237               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
238#endif
239               zdust  = dust(ji,jj) / ( wdust * 30.42 * 0.035 ) * tmask(ji,jj,jk)
240               zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc
241               zscave = zfeequi * zlam1b * zstep
242
243               ! Compute the different ratios for scavenging of iron
244               ! to later allocate scavenged iron to the different organic pools
245               ! ---------------------------------------------------------
246#if  defined key_kriest
247               zdenom1 = xlam1 * trn(ji,jj,jk,jppoc) / zlam1b 
248#else
249               zdenom1 = xlam1 * trn(ji,jj,jk,jppoc) / zlam1b
250               zdenom2 = xlam1 * trn(ji,jj,jk,jpgoc) / zlam1b
251#endif
252
253               !  Increased scavenging for very high iron concentrations found near the coasts
254               !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
255               !  -----------------------------------------------------------
256               zlam1b  = xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
257               zcoag   = zfeequi * zlam1b * zstep
258
259               !  Compute the coagulation of colloidal iron. This parameterization
260               !  could be thought as an equivalent of colloidal pumping.
261               !  It requires certainly some more work as it is very poorly constrained.
262               !  ----------------------------------------------------------------
263               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
264               zlamfac = MIN( 1.  , zlamfac )
265               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
266               zlam1a  = ( 0.369  * 0.3 * trn(ji,jj,jk,jpdoc) + 102.4  * trn(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
267                   &   + ( 114.   * 0.3 * trn(ji,jj,jk,jpdoc) + 5.09E3 * trn(ji,jj,jk,jppoc) )
268#if defined key_kriest
269               zlam1a   = zlam1a + 1E-4 * ( 1. - zlamfac ) * zdep
270               zaggdfea = zlam1a * zstep * zfecoll
271               zaggdfeb = 0.
272               !
273               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
274               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
275#else
276               zlam1b = 3.53E3 *   trn(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) + 1E-4 * ( 1. - zlamfac ) * zdep 
277               zaggdfea = zlam1a * zstep * zfecoll
278               zaggdfeb = zlam1b * zstep * zfecoll
279               !
280               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
281               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
282               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
283#endif
284            END DO
285         END DO
286      END DO
287      !
288      !  Define the bioavailable fraction of iron
289      !  ----------------------------------------
290      IF( ln_fechem ) THEN
291          biron(:,:,:) = MAX( 0., trn(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
292      ELSE
293          biron(:,:,:) = trn(:,:,:,jpfer) 
294      ENDIF
295
296      !  Output of some diagnostics variables
297      !     ---------------------------------
298      IF( ln_diatrc .AND. lk_iomput ) THEN
299         IF( jnt == nrdttrc ) THEN
300            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
301            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
302            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
303            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
304            CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! TL
305            IF( ln_fechem ) THEN
306               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
307               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
308               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
309               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
310            ENDIF
311         ENDIF
312      ENDIF
313
314      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
315         WRITE(charout, FMT="('fechem')")
316         CALL prt_ctl_trc_info(charout)
317         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
318      ENDIF
319      !
320                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
321      IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
322      !
323      IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')
324      !
325   END SUBROUTINE p4z_fechem
326
327
328   SUBROUTINE p4z_fechem_init
329      !!----------------------------------------------------------------------
330      !!                  ***  ROUTINE p4z_fechem_init  ***
331      !!
332      !! ** Purpose :   Initialization of iron chemistry parameters
333      !!
334      !! ** Method  :   Read the nampisfer namelist and check the parameters
335      !!      called at the first timestep
336      !!
337      !! ** input   :   Namelist nampisfer
338      !!
339      !!----------------------------------------------------------------------
340      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand 
341
342      REWIND( numnatp )                     ! read numnatp
343      READ  ( numnatp, nampisfer )
344
345      IF(lwp) THEN                         ! control print
346         WRITE(numout,*) ' '
347         WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
348         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
349         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem =', ln_fechem
350         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar =', ln_ligvar
351         WRITE(numout,*) '    scavenging rate of Iron                   xlam1     =', xlam1
352         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust  =', xlamdust
353         WRITE(numout,*) '    ligand concentration in the ocean         ligand    =', ligand
354      ENDIF
355      !
356      IF( ln_fechem ) THEN
357         ! initialization of some constants used by the complexe chemistry scheme
358         ! ----------------------------------------------------------------------
359         spd = 3600. * 24.
360         con = 1.E9
361         ! LIGAND KINETICS (values from Witter et al. 2000)
362         ! assume Phaeophytin
363         kl1 = 12.2E5  * spd / con
364         kb1 = 12.3E-6 * spd
365         ! Assume DFOB-like for L2
366         kl2 = 19.6e5  * spd / con
367         kb2 = 1.5e-6  * spd
368         ! pcp and remin of Fe3p
369         ks  = 0.075
370         kpr = 0.05
371         ! thermal reduction of Fe3
372         kth = 0.0048 * 24.
373      ENDIF
374      !
375   END SUBROUTINE p4z_fechem_init
376
377#else
378   !!======================================================================
379   !!  Dummy module :                                   No PISCES bio-model
380   !!======================================================================
381CONTAINS
382   SUBROUTINE p4z_fechem                    ! Empty routine
383   END SUBROUTINE p4z_fechem
384#endif 
385
386   !!======================================================================
387END MODULE p4zfechem
Note: See TracBrowser for help on using the repository browser.