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 NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zfechem.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 13.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   !!             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 p4zbc           ! 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_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
28   REAL(wp), PUBLIC ::   xlam1        !: scavenging rate of Iron
29   REAL(wp), PUBLIC ::   xlamdust     !: scavenging rate of Iron by dust
30   REAL(wp), PUBLIC ::   ligand       !: ligand concentration in the ocean
31   REAL(wp), PUBLIC ::   kfep         !: rate constant for nanoparticle formation
32
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
37   !! $Id$
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE p4z_fechem( kt, knt, Kbb, Kmm, Krhs )
43      !!---------------------------------------------------------------------
44      !!                     ***  ROUTINE p4z_fechem  ***
45      !!
46      !! ** Purpose :   Compute remineralization/scavenging of iron
47      !!
48      !! ** Method  :   A simple chemistry model of iron from Aumont and Bopp (2006)
49      !!                based on one ligand and one inorganic form
50      !!---------------------------------------------------------------------
51      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
52      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
53      !
54      INTEGER  ::   ji, jj, jk, jic, jn
55      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
56      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol
57      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
58      REAL(wp) ::   ztrc, zdust
59      REAL(wp) ::   zdenom2
60      REAL(wp) ::   zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2
61      REAL(wp) ::   zrum, zcodel, zargu, zlight
62      REAL(wp) ::   zkox, zkph1, zkph2, zph, zionic, ztligand
63      REAL(wp) ::   za, zb, zc, zkappa1, zkappa2, za0, za1, za2
64      REAL(wp) ::   zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
65      REAL(wp) ::   ztfe, zoxy, zhplus, zxlam
66      REAL(wp) ::   zaggliga, zaggligb
67      REAL(wp) ::   dissol, zligco
68      REAL(wp) :: zrfact2
69      CHARACTER (len=25) :: charout
70      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1
71      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d
72      !!---------------------------------------------------------------------
73      !
74      IF( ln_timing )   CALL timing_start('p4z_fechem')
75      !
76      ! Total ligand concentration : Ligands can be chosen to be constant or variable
77      ! Parameterization from Tagliabue and Voelker (2011)
78      ! -------------------------------------------------
79      IF( ln_ligvar ) THEN
80         ztotlig(:,:,:) =  0.09 * tr(:,:,:,jpdoc,Kbb) * 1E6 + ligand * 1E9
81         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
82      ELSE
83        IF( ln_ligand ) THEN  ;   ztotlig(:,:,:) = tr(:,:,:,jplgw,Kbb) * 1E9
84        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9
85        ENDIF
86      ENDIF
87
88      ! ------------------------------------------------------------
89      !  from Aumont and Bopp (2006)
90      ! This model is based on one ligand and Fe'
91      ! Chemistry is supposed to be fast enough to be at equilibrium
92      ! ------------------------------------------------------------
93      DO_3D_11_11( 1, jpkm1 )
94         zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
95         zkeq            = fekeq(ji,jj,jk)
96         zfesatur        = zTL1(ji,jj,jk) * 1E-9
97         ztfe            = tr(ji,jj,jk,jpfer,Kbb) 
98         ! Fe' is the root of a 2nd order polynom
99         zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
100            &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
101            &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
102         zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
103         zFeL1(ji,jj,jk) = MAX( 0., tr(ji,jj,jk,jpfer,Kbb) * 1E9 - zFe3(ji,jj,jk) )
104      END_3D
105         !
106
107      zdust = 0.         ! if no dust available
108      DO_3D_11_11( 1, jpkm1 )
109         ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water.
110         ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
111         ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
112         ! --------------------------------------------------------------------------------------
113         zhplus  = max( rtrn, hi(ji,jj,jk) )
114         fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
115         &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
116         &         + fesol(ji,jj,jk,5) / zhplus )
117         !
118         zfeequi = zFe3(ji,jj,jk) * 1E-9
119         zhplus  = max( rtrn, hi(ji,jj,jk) )
120         fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
121            &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
122            &         + fesol(ji,jj,jk,5) / zhplus )
123         zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
124         ! precipitation of Fe3+, creation of nanoparticles
125         precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep
126         !
127         ztrc   = ( tr(ji,jj,jk,jppoc,Kbb) + tr(ji,jj,jk,jpgoc,Kbb) + tr(ji,jj,jk,jpcal,Kbb) + tr(ji,jj,jk,jpgsi,Kbb) ) * 1.e6 
128         IF( ll_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) &
129         &  * EXP( -gdept(ji,jj,jk,Kmm) / 540. )
130         IF (ln_ligand) THEN
131            zxlam  = xlam1 * MAX( 1.E-3, EXP(-2 * etot(ji,jj,jk) / 10. ) * (1. - EXP(-2 * tr(ji,jj,jk,jpoxy,Kbb) / 100.E-6 ) ))
132         ELSE
133            zxlam  = xlam1 * 1.0
134         ENDIF
135         zlam1b = 3.e-5 + xlamdust * zdust + zxlam * ztrc
136         zscave = zfeequi * zlam1b * xstep
137
138         ! Compute the different ratios for scavenging of iron
139         ! to later allocate scavenged iron to the different organic pools
140         ! ---------------------------------------------------------
141         zdenom1 = zxlam * tr(ji,jj,jk,jppoc,Kbb) / zlam1b
142         zdenom2 = zxlam * tr(ji,jj,jk,jpgoc,Kbb) / zlam1b
143
144         !  Increased scavenging for very high iron concentrations found near the coasts
145         !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
146         !  -----------------------------------------------------------
147         zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
148         zlamfac = MIN( 1.  , zlamfac )
149         zdep    = MIN( 1., 1000. / gdept(ji,jj,jk,Kmm) )
150         zcoag   = 1E-4 * ( 1. - zlamfac ) * zdep * xstep * tr(ji,jj,jk,jpfer,Kbb)
151
152         !  Compute the coagulation of colloidal iron. This parameterization
153         !  could be thought as an equivalent of colloidal pumping.
154         !  It requires certainly some more work as it is very poorly constrained.
155         !  ----------------------------------------------------------------
156         zlam1a   = ( 0.369  * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) + 102.4  * tr(ji,jj,jk,jppoc,Kbb) ) * xdiss(ji,jj,jk)    &
157             &      + ( 114.   * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) )
158         zaggdfea = zlam1a * xstep * zfecoll
159         !
160         zlam1b   = 3.53E3 * tr(ji,jj,jk,jpgoc,Kbb) * xdiss(ji,jj,jk)
161         zaggdfeb = zlam1b * xstep * zfecoll
162         !
163         tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zscave - zaggdfea - zaggdfeb &
164         &                     - zcoag - precip(ji,jj,jk)
165         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zscave * zdenom1 + zaggdfea
166         tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zscave * zdenom2 + zaggdfeb
167         zscav3d(ji,jj,jk)   = zscave
168         zcoll3d(ji,jj,jk)   = zaggdfea + zaggdfeb
169         !
170      END_3D
171      !
172      !  Define the bioavailable fraction of iron
173      !  ----------------------------------------
174      biron(:,:,:) = tr(:,:,:,jpfer,Kbb) 
175      !
176      IF( ln_ligand ) THEN
177         !
178         DO_3D_11_11( 1, jpkm1 )
179            zlam1a   = ( 0.369  * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) + 102.4  * tr(ji,jj,jk,jppoc,Kbb) ) * xdiss(ji,jj,jk)    &
180                &    + ( 114.   * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) )
181            !
182            zlam1b   = 3.53E3 *   tr(ji,jj,jk,jpgoc,Kbb) * xdiss(ji,jj,jk)
183            zligco   = 0.5 * tr(ji,jj,jk,jplgw,Kmm)
184            zaggliga = zlam1a * xstep * zligco
185            zaggligb = zlam1b * xstep * zligco
186            tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) - zaggliga - zaggligb
187            zlcoll3d(ji,jj,jk)  = zaggliga + zaggligb
188         END_3D
189         !
190         plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( tr(:,:,:,jpfer,Kbb) +rtrn ) ) )
191         !
192      ENDIF
193      !  Output of some diagnostics variables
194      !     ---------------------------------
195      IF( lk_iomput ) THEN
196         IF( knt == nrdttrc ) THEN
197            zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
198            IF( iom_use("Fe3")  )  THEN
199               zFe3(:,:,jpk) = 0.  ;  CALL iom_put("Fe3" , zFe3(:,:,:) * tmask(:,:,:) )   ! Fe3+
200            ENDIF
201            IF( iom_use("FeL1") )  THEN
202              zFeL1(:,:,jpk) = 0.  ;  CALL iom_put("FeL1", zFeL1(:,:,:) * tmask(:,:,:) )   ! FeL1
203            ENDIF
204            IF( iom_use("TL1")  )  THEN
205              zTL1(:,:,jpk) = 0.   ;  CALL iom_put("TL1" , zTL1(:,:,:) * tmask(:,:,:) )   ! TL1
206            ENDIF
207            IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
208            IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:)  * 1e9 * tmask(:,:,:) )   ! biron
209            IF( iom_use("FESCAV") )  THEN
210               zscav3d (:,:,jpk) = 0.  ;  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
211            ENDIF
212            IF( iom_use("FECOLL") ) THEN
213               zcoll3d (:,:,jpk) = 0.  ;   CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
214            ENDIF
215            IF( iom_use("LGWCOLL")) THEN
216               zlcoll3d(:,:,jpk) = 0.  ;  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
217            ENDIF
218          ENDIF
219      ENDIF
220
221      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
222         WRITE(charout, FMT="('fechem')")
223         CALL prt_ctl_trc_info(charout)
224         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
225      ENDIF
226      !
227      IF( ln_timing )   CALL timing_stop('p4z_fechem')
228      !
229   END SUBROUTINE p4z_fechem
230
231
232   SUBROUTINE p4z_fechem_init
233      !!----------------------------------------------------------------------
234      !!                  ***  ROUTINE p4z_fechem_init  ***
235      !!
236      !! ** Purpose :   Initialization of iron chemistry parameters
237      !!
238      !! ** Method  :   Read the nampisfer namelist and check the parameters
239      !!      called at the first timestep
240      !!
241      !! ** input   :   Namelist nampisfer
242      !!
243      !!----------------------------------------------------------------------
244      INTEGER ::   ios   ! Local integer
245      !!
246      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep 
247      !!----------------------------------------------------------------------
248      !
249      IF(lwp) THEN
250         WRITE(numout,*)
251         WRITE(numout,*) 'p4z_rem_init : Initialization of iron chemistry parameters'
252         WRITE(numout,*) '~~~~~~~~~~~~'
253      ENDIF
254      !
255      READ  ( numnatp_ref, nampisfer, IOSTAT = ios, ERR = 901)
256901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisfer in reference namelist' )
257      READ  ( numnatp_cfg, nampisfer, IOSTAT = ios, ERR = 902 )
258902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisfer in configuration namelist' )
259      IF(lwm) WRITE( numonp, nampisfer )
260
261      IF(lwp) THEN                     ! control print
262         WRITE(numout,*) '   Namelist : nampisfer'
263         WRITE(numout,*) '      variable concentration of ligand          ln_ligvar    =', ln_ligvar
264         WRITE(numout,*) '      scavenging rate of Iron                   xlam1        =', xlam1
265         WRITE(numout,*) '      scavenging rate of Iron by dust           xlamdust     =', xlamdust
266         WRITE(numout,*) '      ligand concentration in the ocean         ligand       =', ligand
267         WRITE(numout,*) '      rate constant for nanoparticle formation  kfep         =', kfep
268      ENDIF
269      !
270   END SUBROUTINE p4z_fechem_init
271   
272   !!======================================================================
273END MODULE p4zfechem
Note: See TracBrowser for help on using the repository browser.