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/OCE/ISF – 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.

Location:
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfcavmlt.F90

    r12077 r12340  
    3131   PUBLIC   isfcav_mlt 
    3232 
     33   !! * Substitutions 
     34#  include "do_loop_substitute.h90" 
    3335   !!---------------------------------------------------------------------- 
    3436   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    208210      ! compute upward heat flux zhtflx and upward water flux zwflx 
    209211      ! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24) 
    210       DO jj = 1, jpj 
    211          DO ji = 1, jpi 
    212             ! 
    213             ! compute coeficient to solve the 2nd order equation 
    214             zeps1 = rau0_rcp * pgt(ji,jj) 
    215             zeps2 = rLfusisf * rau0 * pgs(ji,jj) 
    216             zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps) 
    217             zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj) 
    218             zeps6 = zeps4 - pttbl(ji,jj) 
    219             zeps7 = zeps4 - rtsurf 
    220             ! 
    221             ! solve the 2nd order equation to find zsfrz 
    222             zaqe  = risf_lamb1 * (zeps1 + zeps3) 
    223             zaqer = 0.5_wp / MIN(zaqe,-zeps) 
    224             zbqe  = zeps1 * zeps6 + zeps3 * zeps7 - zeps2 
    225             zcqe  = zeps2 * pstbl(ji,jj) 
    226             zdis  = zbqe * zbqe - 4.0_wp * zaqe * zcqe                
    227             ! 
    228             ! Presumably zdis can never be negative because gammas is very small compared to gammat 
    229             zsfrz=(-zbqe - SQRT(zdis)) * zaqer 
    230             IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe + SQRT(zdis)) * zaqer  ! check this if this if is needed 
    231             ! 
    232             ! compute t freeze (eq. 25) 
    233             ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz 
    234             ! 
    235             ! thermal driving 
    236             zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) ) 
    237             ! 
    238             ! compute the upward water and heat flux (eq. 24 and eq. 26) 
    239             pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out) 
    240             pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out) 
    241             pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out) 
    242             ! 
    243             zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf ) 
    244             ! 
    245          END DO 
    246       END DO 
     212      DO_2D_11_11 
     213         ! 
     214         ! compute coeficient to solve the 2nd order equation 
     215         zeps1 = rau0_rcp * pgt(ji,jj) 
     216         zeps2 = rLfusisf * rau0 * pgs(ji,jj) 
     217         zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps) 
     218         zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj) 
     219         zeps6 = zeps4 - pttbl(ji,jj) 
     220         zeps7 = zeps4 - rtsurf 
     221         ! 
     222         ! solve the 2nd order equation to find zsfrz 
     223         zaqe  = risf_lamb1 * (zeps1 + zeps3) 
     224         zaqer = 0.5_wp / MIN(zaqe,-zeps) 
     225         zbqe  = zeps1 * zeps6 + zeps3 * zeps7 - zeps2 
     226         zcqe  = zeps2 * pstbl(ji,jj) 
     227         zdis  = zbqe * zbqe - 4.0_wp * zaqe * zcqe                
     228         ! 
     229         ! Presumably zdis can never be negative because gammas is very small compared to gammat 
     230         zsfrz=(-zbqe - SQRT(zdis)) * zaqer 
     231         IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe + SQRT(zdis)) * zaqer  ! check this if this if is needed 
     232         ! 
     233         ! compute t freeze (eq. 25) 
     234         ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz 
     235         ! 
     236         ! thermal driving 
     237         zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) ) 
     238         ! 
     239         ! compute the upward water and heat flux (eq. 24 and eq. 26) 
     240         pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out) 
     241         pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out) 
     242         pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out) 
     243         ! 
     244         zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf ) 
     245         ! 
     246      END_2D 
    247247      ! 
    248248      ! output conductive heat flux through the ice 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfcpl.F90

    r12077 r12340  
    4040   END TYPE 
    4141   ! 
     42   !! * Substitutions 
     43#  include "do_loop_substitute.h90" 
    4244   !!---------------------------------------------------------------------- 
    4345   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    360362      ! ----------------------------------------------------------------------------------------- 
    361363      ! case we open a cell but no neigbour cells available to get an estimate of T and S 
    362       DO jk = 1,jpk-1 
    363          DO jj = 1,jpj 
    364             DO ji = 1,jpi 
    365                IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp)              & 
    366                   &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     & 
    367                   &                          try increase nn_drown or activate XXXX & 
    368                   &                         in your domain cfg computation'         ) 
    369             END DO 
    370          END DO 
    371       END DO 
     364      DO_3D_11_11( 1,jpk-1 ) 
     365         IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp)              & 
     366            &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     & 
     367            &                          try increase nn_drown or activate XXXX & 
     368            &                         in your domain cfg computation'         ) 
     369      END_3D 
    372370      !  
    373371   END SUBROUTINE isfcpl_tra 
    … …  
    404402      DO jk = 1, jpk                                 ! Horizontal slab 
    405403         ! 1.1: get volume flux before coupling (>0 out) 
    406          DO jj = 2, jpjm1 
    407             DO ji = 2, jpim1 
    408                zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * uu(ji-1,jj  ,jk,Kmm)    & 
    409                   &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vv(ji  ,jj-1,jk,Kmm)  ) & 
    410                   &                * ztmask_b(ji,jj,jk) 
    411             END DO 
    412          ENDDO 
     404         DO_2D_00_00 
     405            zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * uu(ji-1,jj  ,jk,Kmm)    & 
     406               &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vv(ji  ,jj-1,jk,Kmm)  ) & 
     407               &                * ztmask_b(ji,jj,jk) 
     408         END_2D 
    413409         ! 
    414410         ! 1.2: get volume flux after coupling (>0 out) 
    … …  
    418414         vv(:,:,jk,Kmm) = vv(:,:,jk,Kmm) * vmask(:,:,jk) 
    419415         ! compute volume flux divergence after coupling 
    420          DO jj = 2, jpjm1 
    421             DO ji = 2, jpim1 
    422                zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * uu(ji-1,jj  ,jk,Kmm)    & 
    423                   &                 + e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * vv(ji  ,jj-1,jk,Kmm)  ) & 
    424                   &               * tmask(ji,jj,jk) 
    425             END DO 
    426          ENDDO 
     416         DO_2D_00_00 
     417            zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj  ) * e3u(ji-1,jj  ,jk,Kmm) * uu(ji-1,jj  ,jk,Kmm)    & 
     418               &                 + e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) - e1v(ji  ,jj-1) * e3v(ji  ,jj-1,jk,Kmm) * vv(ji  ,jj-1,jk,Kmm)  ) & 
     419               &               * tmask(ji,jj,jk) 
     420         END_2D 
    427421         ! 
    428422         ! 1.3: get 3d volume flux difference (before - after cpl) (>0 out) 
    … …  
    433427      ! 2.0: include the contribution of the vertical velocity in the volume flux correction 
    434428      ! 
    435       DO jj = 2, jpjm1 
    436          DO ji = 2, jpim1 
    437             ! 
    438             ikt = mikt(ji,jj) 
    439             IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN 
    440                risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1))  ! test sign 
    441             ENDIF 
    442             ! 
    443          END DO 
    444       ENDDO 
     429      DO_2D_00_00 
     430         ! 
     431         ikt = mikt(ji,jj) 
     432         IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN 
     433            risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1))  ! test sign 
     434         ENDIF 
     435         ! 
     436      END_2D 
    445437      ! 
    446438      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. ) 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfdiags.F90

    r12077 r12340  
    2424   PUBLIC   isf_diags_flx 
    2525 
     26   !! * Substitutions 
     27#  include "do_loop_substitute.h90" 
    2628   !!---------------------------------------------------------------------- 
    2729   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    98100      zvar3d(:,:,:) = 0._wp 
    99101      ! 
    100       DO jj = 1,jpj 
    101          DO ji = 1,jpi 
    102             ikt = ktop(ji,jj) 
    103             ikb = kbot(ji,jj) 
    104             DO jk = ikt, ikb - 1 
    105                zvar3d(ji,jj,jk) = zvar2d(ji,jj) * e3t(ji,jj,jk,Kmm) 
    106             END DO 
    107             zvar3d(ji,jj,ikb) = zvar2d(ji,jj) * e3t(ji,jj,ikb,Kmm) * pfrac(ji,jj) 
     102      DO_2D_11_11 
     103         ikt = ktop(ji,jj) 
     104         ikb = kbot(ji,jj) 
     105         DO jk = ikt, ikb - 1 
     106            zvar3d(ji,jj,jk) = zvar2d(ji,jj) * e3t(ji,jj,jk,Kmm) 
    108107         END DO 
    109       END DO 
     108         zvar3d(ji,jj,ikb) = zvar2d(ji,jj) * e3t(ji,jj,ikb,Kmm) * pfrac(ji,jj) 
     109      END_2D 
    110110      ! 
    111111      CALL iom_put( TRIM(cdvar) , zvar3d(:,:,:)) 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfdynatf.F90

    r12077 r12340  
    2323 
    2424   PUBLIC isf_dynatf 
     25   !! * Substitutions 
     26#  include "do_loop_substitute.h90" 
    2527 
    2628CONTAINS 
    … …  
    7880      ! 
    7981      ! add the increment in the tbl 
    80       DO jk = 1, jpkm1 
    81          DO jj = 1, jpj 
    82             DO ji = 1, jpi 
    83                IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN 
    84                   pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm) 
    85                ELSEIF ( jk == kbot(ji,jj) ) THEN 
    86                   pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm) * pfrac(ji,jj) 
    87                ENDIF 
    88             END DO 
    89          END DO 
    90       END DO 
     82      DO_3D_11_11( 1, jpkm1 ) 
     83         IF( ktop(ji,jj) <= jk .AND. jk < kbot(ji,jj)  ) THEN 
     84            pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm) 
     85         ELSEIF ( jk == kbot(ji,jj) ) THEN 
     86            pe3t_f(ji,jj,jk) = pe3t_f(ji,jj,jk) - zfwfinc(ji,jj) * e3t(ji,jj,jk,Kmm) * pfrac(ji,jj) 
     87         ENDIF 
     88      END_3D 
    9189      ! 
    9290   END SUBROUTINE isf_dynatf_mlt 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfhdiv.F90

    r12077 r12340  
    2424 
    2525   PUBLIC isf_hdiv 
     26   !! * Substitutions 
     27#  include "do_loop_substitute.h90" 
    2628 
    2729CONTAINS 
    … …  
    9799      ! 
    98100      ! update divergence at each level affected by ice shelf top boundary layer 
    99       DO jj = 1,jpj 
    100          DO ji = 1,jpi 
    101             ikt = ktop(ji,jj) 
    102             ikb = kbot(ji,jj) 
    103             ! level fully include in the ice shelf boundary layer 
    104             DO jk = ikt, ikb - 1 
    105                phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zhdiv(ji,jj) 
    106             END DO 
    107             ! level partially include in ice shelf boundary layer  
    108             phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zhdiv(ji,jj) * pfrac(ji,jj) 
     101      DO_2D_11_11 
     102         ikt = ktop(ji,jj) 
     103         ikb = kbot(ji,jj) 
     104         ! level fully include in the ice shelf boundary layer 
     105         DO jk = ikt, ikb - 1 
     106            phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zhdiv(ji,jj) 
    109107         END DO 
    110       END DO 
     108         ! level partially include in ice shelf boundary layer  
     109         phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zhdiv(ji,jj) * pfrac(ji,jj) 
     110      END_2D 
    111111      ! 
    112112   END SUBROUTINE isf_hdiv_mlt 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfload.F90

    r12077 r12340  
    2424 
    2525   PUBLIC isf_load 
     26   !! * Substitutions 
     27#  include "do_loop_substitute.h90" 
    2628 
    2729CONTAINS 
    … …  
    9193      !                                !- Surface value + ice shelf gradient 
    9294      pisfload(:,:) = 0._wp                       ! compute pressure due to ice shelf load  
    93       DO jj = 1, jpj                         ! (used to compute hpgi/j for all the level from 1 to miku/v) 
    94          DO ji = 1, jpi                      ! divided by 2 later 
    95             ikt = mikt(ji,jj) 
     95      DO_2D_11_11 
     96         ikt = mikt(ji,jj) 
     97         ! 
     98         IF ( ikt > 1 ) THEN 
    9699            ! 
    97             IF ( ikt > 1 ) THEN 
    98                ! 
    99                ! top layer of the ice shelf 
    100                pisfload(ji,jj) = pisfload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) 
    101                ! 
    102                ! core layers of the ice shelf 
    103                DO jk = 2, ikt-1 
    104                   pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) 
    105                END DO 
    106                ! 
    107                ! deepest part of the ice shelf (between deepest T point and ice/ocean interface 
    108                pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
    109                   &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 
    110                ! 
    111             END IF 
    112          END DO 
    113       END DO 
     100            ! top layer of the ice shelf 
     101            pisfload(ji,jj) = pisfload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) 
     102            ! 
     103            ! core layers of the ice shelf 
     104            DO jk = 2, ikt-1 
     105               pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) 
     106            END DO 
     107            ! 
     108            ! deepest part of the ice shelf (between deepest T point and ice/ocean interface 
     109            pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
     110               &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 
     111            ! 
     112         END IF 
     113      END_2D 
    114114      ! 
    115115   END SUBROUTINE isf_load_uniform 

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isftbl.F90

    r12336 r12340  
    2323 
    2424   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot 
     25   !! * Substitutions 
     26#  include "do_loop_substitute.h90" 
    2527 
    2628CONTAINS 
    … …  
    7072         ! compute tbl property at T point 
    7173         pvarout(1,:) = 0._wp 
    72          DO jj = 1, jpj 
    73             DO ji = 2, jpi 
    74                pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj)) 
    75             END DO 
    76          END DO 
     74         DO_2D_11_01 
     75            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj)) 
     76         END_2D 
    7777         ! lbclnk not needed as a final communication is done after the computation of fwf 
    7878         !  
    … …  
    9090         ! pvarout is an averaging of wet point 
    9191         pvarout(:,1) = 0._wp 
    92          DO jj = 2, jpj 
    93             DO ji = 1, jpi 
    94                pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1)) 
    95             END DO 
    96          END DO 
     92         DO_2D_01_11 
     93            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1)) 
     94         END_2D 
    9795         ! lbclnk not needed as a final communication is done after the computation of fwf 
    9896         ! 
    … …  
    128126      ! 
    129127      ! compute tbl top.bottom level and thickness 
    130       DO jj = 1,jpj 
    131          DO ji = 1,jpi 
    132             ! 
    133             ! tbl top/bottom indices initialisation 
    134             ikt = ktop(ji,jj) ; ikb = kbot(ji,jj) 
    135             ! 
    136             ! level fully include in the ice shelf boundary layer 
    137             pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj) 
    138             ! 
    139             ! level partially include in ice shelf boundary layer  
    140             pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj) 
    141             ! 
    142          END DO 
    143       END DO 
     128      DO_2D_11_11 
     129         ! 
     130         ! tbl top/bottom indices initialisation 
     131         ikt = ktop(ji,jj) ; ikb = kbot(ji,jj) 
     132         ! 
     133         ! level fully include in the ice shelf boundary layer 
     134         pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj) 
     135         ! 
     136         ! level partially include in ice shelf boundary layer  
     137         pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj) 
     138         ! 
     139      END_2D 
    144140 
    145141   END SUBROUTINE isf_tbl_avg 
    … …  
    168164      ! 
    169165      ! get htbl 
    170       DO jj = 1,jpj 
    171          DO ji = 1,jpi 
    172             ! 
    173             ! tbl top/bottom indices initialisation 
    174             ikt = ktop(ji,jj) 
    175             ! 
    176             ! limit the tbl to water thickness. 
    177             phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) ) 
    178             ! 
    179             ! thickness of boundary layer must be at least the top level thickness 
    180             phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) ) 
    181             ! 
    182          END DO 
    183       END DO 
     166      DO_2D_11_11 
     167         ! 
     168         ! tbl top/bottom indices initialisation 
     169         ikt = ktop(ji,jj) 
     170         ! 
     171         ! limit the tbl to water thickness. 
     172         phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) ) 
     173         ! 
     174         ! thickness of boundary layer must be at least the top level thickness 
     175         phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) ) 
     176         ! 
     177      END_2D 
    184178      ! 
    185179      ! get ktbl 
    … …  
    187181      ! 
    188182      ! get pfrac 
    189       DO jj = 1,jpj 
    190          DO ji = 1,jpi 
    191             ! 
    192             ! tbl top/bottom indices initialisation 
    193             ikt = ktop(ji,jj) ; ikb = kbot(ji,jj) 
    194             ! 
    195             ! proportion of the bottom cell included in ice shelf boundary layer  
    196             pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb) 
    197             ! 
    198          END DO 
    199       END DO 
     183      DO_2D_11_11 
     184         ! 
     185         ! tbl top/bottom indices initialisation 
     186         ikt = ktop(ji,jj) ; ikb = kbot(ji,jj) 
     187         ! 
     188         ! proportion of the bottom cell included in ice shelf boundary layer  
     189         pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb) 
     190         ! 
     191      END_2D 
    200192      ! 
    201193   END SUBROUTINE isf_tbl_lvl 
    … …  
    223215      ! 
    224216      ! get ktbl 
    225       DO jj = 1,jpj 
    226          DO ji = 1,jpi 
    227             ! 
    228             ! determine the deepest level influenced by the boundary layer 
    229             ikt = ktop(ji,jj) 
    230             ikb = ikt 
    231             DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ;  ikb = ikb + 1 ;  END DO 
    232             kbot(ji,jj) = ikb - 1 
    233             ! 
    234          END DO 
    235       END DO 
     217      DO_2D_11_11 
     218         ! 
     219         ! determine the deepest level influenced by the boundary layer 
     220         ikt = ktop(ji,jj) 
     221         ikb = ikt 
     222         DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ;  ikb = ikb + 1 ;  END DO 
     223         kbot(ji,jj) = ikb - 1 
     224         ! 
     225      END_2D 
    236226      ! 
    237227   END SUBROUTINE isf_tbl_kbot 
    … …  
    259249      ! test: this routine run with pdep = 0 should return 1 
    260250      ! 
    261       DO jj = 1, jpj 
    262          DO ji = 1, jpi 
    263             ! comput ktop 
    264             ikt = 2 
    265             DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ;  ikt = ikt + 1 ;  END DO 
    266             ktop(ji,jj) = ikt - 1 
    267             ! 
    268             ! update pdep 
    269             pdep(ji,jj) = gdepw_0(ji,jj,ktop(ji,jj)) 
    270          END DO 
    271       END DO 
     251      DO_2D_11_11 
     252         ! comput ktop 
     253         ikt = 2 
     254         DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ;  ikt = ikt + 1 ;  END DO 
     255         ktop(ji,jj) = ikt - 1 
     256         ! 
     257         ! update pdep 
     258         pdep(ji,jj) = gdepw_0(ji,jj,ktop(ji,jj)) 
     259      END_2D 
    272260      ! 
    273261   END SUBROUTINE isf_tbl_ktop 
Note: See TracChangeset for help on using the changeset viewer.