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.
Changeset 12340 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2020-01-27T15:31:53+01:00 (4 years ago)
Author:
acc
Message:

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zsed.F90

    r12258 r12340  
    3737   REAL(wp), SAVE :: sedsilfrac, sedcalfrac 
    3838 
     39   !! * Substitutions 
     40#  include "do_loop_substitute.h90" 
    3941   !!---------------------------------------------------------------------- 
    4042   !! NEMO/TOP 4.0 , NEMO Consortium (2018) 
    … …  
    9193         ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 
    9294         ! -------------------------------------------------------------------- 
    93          DO jj = 1, jpj 
    94             DO ji = 1, jpi 
    95                ikt  = mbkt(ji,jj) 
    96                zdep = e3t(ji,jj,ikt,Kmm) / xstep 
    97                zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 
    98                zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 
    99             END DO 
    100          END DO 
     95         DO_2D_11_11 
     96            ikt  = mbkt(ji,jj) 
     97            zdep = e3t(ji,jj,ikt,Kmm) / xstep 
     98            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 
     99            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 
     100         END_2D 
    101101 
    102102         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 
    103103         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 
    104104         ! ------------------------------------------------------- 
    105          DO jj = 1, jpj 
    106             DO ji = 1, jpi 
    107               IF( tmask(ji,jj,1) == 1 ) THEN 
    108                  ikt = mbkt(ji,jj) 
    109                  zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
    110                    &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
    111                  zflx  = LOG10( MAX( 1E-3, zflx ) ) 
    112                  zo2   = LOG10( MAX( 10. , tr(ji,jj,ikt,jpoxy,Kbb) * 1E6 ) ) 
    113                  zno3  = LOG10( MAX( 1.  , tr(ji,jj,ikt,jpno3,Kbb) * 1E6 * rno3 ) ) 
    114                  zdep  = LOG10( gdepw(ji,jj,ikt+1,Kmm) ) 
    115                  zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
    116                    &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
    117                  zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
    118                    ! 
    119                  zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
    120                    &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6 
    121                  zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
    122               ENDIF 
    123             END DO 
    124          END DO  
     105         DO_2D_11_11 
     106           IF( tmask(ji,jj,1) == 1 ) THEN 
     107              ikt = mbkt(ji,jj) 
     108              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
     109                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
     110              zflx  = LOG10( MAX( 1E-3, zflx ) ) 
     111              zo2   = LOG10( MAX( 10. , tr(ji,jj,ikt,jpoxy,Kbb) * 1E6 ) ) 
     112              zno3  = LOG10( MAX( 1.  , tr(ji,jj,ikt,jpno3,Kbb) * 1E6 * rno3 ) ) 
     113              zdep  = LOG10( gdepw(ji,jj,ikt+1,Kmm) ) 
     114              zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
     115                &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
     116              zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
     117                ! 
     118              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
     119                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6 
     120              zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
     121           ENDIF 
     122         END_2D 
    125123         ! 
    126124      ENDIF 
    … …  
    131129      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac 
    132130 
    133       DO jj = 1, jpj 
    134          DO ji = 1, jpi 
     131      DO_2D_11_11 
     132         ikt  = mbkt(ji,jj) 
     133         zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     134         zwsc = zwsbio4(ji,jj) * zdep 
     135         zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc 
     136         zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc 
     137         ! 
     138         tr(ji,jj,ikt,jpgsi,Krhs) = tr(ji,jj,ikt,jpgsi,Krhs) - zsiloss 
     139         tr(ji,jj,ikt,jpcal,Krhs) = tr(ji,jj,ikt,jpcal,Krhs) - zcaloss 
     140      END_2D 
     141      ! 
     142      IF( .NOT.lk_sed ) THEN 
     143         DO_2D_11_11 
    135144            ikt  = mbkt(ji,jj) 
    136145            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    … …  
    138147            zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc 
    139148            zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc 
     149            tr(ji,jj,ikt,jpsil,Krhs) = tr(ji,jj,ikt,jpsil,Krhs) + zsiloss * zrivsil  
    140150            ! 
    141             tr(ji,jj,ikt,jpgsi,Krhs) = tr(ji,jj,ikt,jpgsi,Krhs) - zsiloss 
    142             tr(ji,jj,ikt,jpcal,Krhs) = tr(ji,jj,ikt,jpcal,Krhs) - zcaloss 
    143          END DO 
    144       END DO 
    145       ! 
    146       IF( .NOT.lk_sed ) THEN 
    147          DO jj = 1, jpj 
    148             DO ji = 1, jpi 
    149                ikt  = mbkt(ji,jj) 
    150                zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    151                zwsc = zwsbio4(ji,jj) * zdep 
    152                zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc 
    153                zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc 
    154                tr(ji,jj,ikt,jpsil,Krhs) = tr(ji,jj,ikt,jpsil,Krhs) + zsiloss * zrivsil  
    155                ! 
    156                zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
    157                zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    158                zrivalk  = sedcalfrac * zfactcal 
    159                tr(ji,jj,ikt,jptal,Krhs) =  tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk * 2.0 
    160                tr(ji,jj,ikt,jpdic,Krhs) =  tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk 
    161                zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t(ji,jj,ikt,Kmm)  
    162                zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t(ji,jj,ikt,Kmm)  
    163             END DO 
    164          END DO 
    165       ENDIF 
    166       ! 
    167       DO jj = 1, jpj 
    168          DO ji = 1, jpi 
     151            zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
     152            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     153            zrivalk  = sedcalfrac * zfactcal 
     154            tr(ji,jj,ikt,jptal,Krhs) =  tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk * 2.0 
     155            tr(ji,jj,ikt,jpdic,Krhs) =  tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk 
     156            zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t(ji,jj,ikt,Kmm)  
     157            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t(ji,jj,ikt,Kmm)  
     158         END_2D 
     159      ENDIF 
     160      ! 
     161      DO_2D_11_11 
     162         ikt  = mbkt(ji,jj) 
     163         zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     164         zws4 = zwsbio4(ji,jj) * zdep 
     165         zws3 = zwsbio3(ji,jj) * zdep 
     166         tr(ji,jj,ikt,jpgoc,Krhs) = tr(ji,jj,ikt,jpgoc,Krhs) - tr(ji,jj,ikt,jpgoc,Kbb) * zws4  
     167         tr(ji,jj,ikt,jppoc,Krhs) = tr(ji,jj,ikt,jppoc,Krhs) - tr(ji,jj,ikt,jppoc,Kbb) * zws3 
     168         tr(ji,jj,ikt,jpbfe,Krhs) = tr(ji,jj,ikt,jpbfe,Krhs) - tr(ji,jj,ikt,jpbfe,Kbb) * zws4 
     169         tr(ji,jj,ikt,jpsfe,Krhs) = tr(ji,jj,ikt,jpsfe,Krhs) - tr(ji,jj,ikt,jpsfe,Kbb) * zws3 
     170      END_2D 
     171      ! 
     172      IF( ln_p5z ) THEN 
     173         DO_2D_11_11 
    169174            ikt  = mbkt(ji,jj) 
    170175            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    171176            zws4 = zwsbio4(ji,jj) * zdep 
    172177            zws3 = zwsbio3(ji,jj) * zdep 
    173             tr(ji,jj,ikt,jpgoc,Krhs) = tr(ji,jj,ikt,jpgoc,Krhs) - tr(ji,jj,ikt,jpgoc,Kbb) * zws4  
    174             tr(ji,jj,ikt,jppoc,Krhs) = tr(ji,jj,ikt,jppoc,Krhs) - tr(ji,jj,ikt,jppoc,Kbb) * zws3 
    175             tr(ji,jj,ikt,jpbfe,Krhs) = tr(ji,jj,ikt,jpbfe,Krhs) - tr(ji,jj,ikt,jpbfe,Kbb) * zws4 
    176             tr(ji,jj,ikt,jpsfe,Krhs) = tr(ji,jj,ikt,jpsfe,Krhs) - tr(ji,jj,ikt,jpsfe,Kbb) * zws3 
    177          END DO 
    178       END DO 
    179       ! 
    180       IF( ln_p5z ) THEN 
    181          DO jj = 1, jpj 
    182             DO ji = 1, jpi 
    183                ikt  = mbkt(ji,jj) 
    184                zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    185                zws4 = zwsbio4(ji,jj) * zdep 
    186                zws3 = zwsbio3(ji,jj) * zdep 
    187                tr(ji,jj,ikt,jpgon,Krhs) = tr(ji,jj,ikt,jpgon,Krhs) - tr(ji,jj,ikt,jpgon,Kbb) * zws4 
    188                tr(ji,jj,ikt,jppon,Krhs) = tr(ji,jj,ikt,jppon,Krhs) - tr(ji,jj,ikt,jppon,Kbb) * zws3 
    189                tr(ji,jj,ikt,jpgop,Krhs) = tr(ji,jj,ikt,jpgop,Krhs) - tr(ji,jj,ikt,jpgop,Kbb) * zws4 
    190                tr(ji,jj,ikt,jppop,Krhs) = tr(ji,jj,ikt,jppop,Krhs) - tr(ji,jj,ikt,jppop,Kbb) * zws3 
    191             END DO 
    192          END DO 
     178            tr(ji,jj,ikt,jpgon,Krhs) = tr(ji,jj,ikt,jpgon,Krhs) - tr(ji,jj,ikt,jpgon,Kbb) * zws4 
     179            tr(ji,jj,ikt,jppon,Krhs) = tr(ji,jj,ikt,jppon,Krhs) - tr(ji,jj,ikt,jppon,Kbb) * zws3 
     180            tr(ji,jj,ikt,jpgop,Krhs) = tr(ji,jj,ikt,jpgop,Krhs) - tr(ji,jj,ikt,jpgop,Kbb) * zws4 
     181            tr(ji,jj,ikt,jppop,Krhs) = tr(ji,jj,ikt,jppop,Krhs) - tr(ji,jj,ikt,jppop,Kbb) * zws3 
     182         END_2D 
    193183      ENDIF 
    194184 
    … …  
    196186         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after 
    197187         ! denitrification in the sediments. Not very clever, but simpliest option. 
    198          DO jj = 1, jpj 
    199             DO ji = 1, jpi 
    200                ikt  = mbkt(ji,jj) 
    201                zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    202                zws4 = zwsbio4(ji,jj) * zdep 
    203                zws3 = zwsbio3(ji,jj) * zdep 
    204                zrivno3 = 1. - zbureff(ji,jj) 
    205                zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3 
    206                zpdenit  = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
    207                z1pdenit = zwstpoc * zrivno3 - zpdenit 
    208                zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
    209                tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit 
    210                tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit 
    211                tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit 
    212                tr(ji,jj,ikt,jpno3,Krhs) = tr(ji,jj,ikt,jpno3,Krhs) - rdenit * zpdenit 
    213                tr(ji,jj,ikt,jpoxy,Krhs) = tr(ji,jj,ikt,jpoxy,Krhs) - zolimit * o2ut 
    214                tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + rno3 * (zolimit + (1.+rdenit) * zpdenit ) 
    215                tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zpdenit + zolimit  
    216                sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm) 
    217                zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm) 
    218                IF( ln_p5z ) THEN 
    219                   zwstpop              = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3 
    220                   zwstpon              = tr(ji,jj,ikt,jpgon,Kbb) * zws4 + tr(ji,jj,ikt,jppon,Kbb) * zws3 
    221                   tr(ji,jj,ikt,jpdon,Krhs) = tr(ji,jj,ikt,jpdon,Krhs) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn) 
    222                   tr(ji,jj,ikt,jpdop,Krhs) = tr(ji,jj,ikt,jpdop,Krhs) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn) 
    223                ENDIF 
    224             END DO 
    225          END DO 
     188         DO_2D_11_11 
     189            ikt  = mbkt(ji,jj) 
     190            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     191            zws4 = zwsbio4(ji,jj) * zdep 
     192            zws3 = zwsbio3(ji,jj) * zdep 
     193            zrivno3 = 1. - zbureff(ji,jj) 
     194            zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3 
     195            zpdenit  = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
     196            z1pdenit = zwstpoc * zrivno3 - zpdenit 
     197            zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
     198            tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit 
     199            tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit 
     200            tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit 
     201            tr(ji,jj,ikt,jpno3,Krhs) = tr(ji,jj,ikt,jpno3,Krhs) - rdenit * zpdenit 
     202            tr(ji,jj,ikt,jpoxy,Krhs) = tr(ji,jj,ikt,jpoxy,Krhs) - zolimit * o2ut 
     203            tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + rno3 * (zolimit + (1.+rdenit) * zpdenit ) 
     204            tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zpdenit + zolimit  
     205            sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm) 
     206            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm) 
     207            IF( ln_p5z ) THEN 
     208               zwstpop              = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3 
     209               zwstpon              = tr(ji,jj,ikt,jpgon,Kbb) * zws4 + tr(ji,jj,ikt,jppon,Kbb) * zws3 
     210               tr(ji,jj,ikt,jpdon,Krhs) = tr(ji,jj,ikt,jpdon,Krhs) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn) 
     211               tr(ji,jj,ikt,jpdop,Krhs) = tr(ji,jj,ikt,jpdop,Krhs) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn) 
     212            ENDIF 
     213         END_2D 
    226214       ENDIF 
    227215 
    … …  
    235223      ENDDO 
    236224      IF( ln_p4z ) THEN 
    237          DO jk = 1, jpkm1 
    238             DO jj = 1, jpj 
    239                DO ji = 1, jpi 
    240                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    241                   ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
    242                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    243                   !       Potential nitrogen fixation dependant on temperature and iron 
    244                   xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
    245                   xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
    246                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    247                   IF( zlim <= 0.1 )   zlim = 0.01 
    248                   zfact = zlim * rfact2 
    249                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    250                   ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
    251                   ztrdp = ztrpo4(ji,jj,jk) 
    252                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    253                END DO 
    254             END DO 
    255          END DO 
     225         DO_3D_11_11( 1, jpkm1 ) 
     226            !                      ! Potential nitrogen fixation dependant on temperature and iron 
     227            ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
     228            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
     229            !       Potential nitrogen fixation dependant on temperature and iron 
     230            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
     231            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
     232            zlim = ( 1.- xdiano3 - xdianh4 ) 
     233            IF( zlim <= 0.1 )   zlim = 0.01 
     234            zfact = zlim * rfact2 
     235            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
     236            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
     237            ztrdp = ztrpo4(ji,jj,jk) 
     238            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
     239         END_3D 
    256240      ELSE       ! p5z 
    257          DO jk = 1, jpkm1 
    258             DO jj = 1, jpj 
    259                DO ji = 1, jpi 
    260                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    261                   ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
    262                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    263                   !       Potential nitrogen fixation dependant on temperature and iron 
    264                   xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
    265                   xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
    266                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    267                   IF( zlim <= 0.1 )   zlim = 0.01 
    268                   zfact = zlim * rfact2 
    269                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    270                   ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
    271                   ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk)) 
    272                   ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
    273                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    274                END DO 
    275             END DO 
    276          END DO 
     241         DO_3D_11_11( 1, jpkm1 ) 
     242            !                      ! Potential nitrogen fixation dependant on temperature and iron 
     243            ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
     244            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
     245            !       Potential nitrogen fixation dependant on temperature and iron 
     246            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
     247            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
     248            zlim = ( 1.- xdiano3 - xdianh4 ) 
     249            IF( zlim <= 0.1 )   zlim = 0.01 
     250            zfact = zlim * rfact2 
     251            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
     252            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
     253            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk)) 
     254            ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
     255            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
     256         END_3D 
    277257      ENDIF 
    278258 
    … …  
    280260      ! ---------------------------------------- 
    281261      IF( ln_p4z ) THEN 
    282          DO jk = 1, jpkm1 
    283             DO jj = 1, jpj 
    284                DO ji = 1, jpi 
    285                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    286                   tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
    287                   tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
    288                   tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zfact * 2.0 / 3.0 
    289                   tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
    290                   tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    291                   tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    292                   tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    293                   tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 
    294                   tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    295                   tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    296                   tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    297                   tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) & 
    298                   &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep 
    299               END DO 
    300             END DO  
    301          END DO 
     262         DO_3D_11_11( 1, jpkm1 ) 
     263            zfact = nitrpot(ji,jj,jk) * nitrfix 
     264            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
     265            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
     266            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zfact * 2.0 / 3.0 
     267            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
     268            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
     269            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
     270            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
     271            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 
     272            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
     273            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
     274            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     275            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) & 
     276            &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep 
     277         END_3D 
    302278      ELSE    ! p5z 
    303          DO jk = 1, jpkm1 
    304             DO jj = 1, jpj 
    305                DO ji = 1, jpi 
    306                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    307                   tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
    308                   tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
    309                   tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
    310                   &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    311                   tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0 
    312                   tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
    313                   tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  & 
    314                   &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
    315                   &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    316                   tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    317                   tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0 
    318                   tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
    319                   tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    320                   tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) + zfact * 1.0 / 3.0 * 1.0 /3.0 
    321                   tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
    322                   tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    323                   tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0  
    324                   tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    325                   tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    326                   tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    327               END DO 
    328             END DO  
    329          END DO 
     279         DO_3D_11_11( 1, jpkm1 ) 
     280            zfact = nitrpot(ji,jj,jk) * nitrfix 
     281            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
     282            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
     283            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
     284            &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     285            tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0 
     286            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
     287            tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  & 
     288            &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
     289            &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     290            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
     291            tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0 
     292            tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
     293            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
     294            tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) + zfact * 1.0 / 3.0 * 1.0 /3.0 
     295            tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
     296            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
     297            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0  
     298            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
     299            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
     300            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     301         END_3D 
    330302         ! 
    331303      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.