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

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

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