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/ICE/icedyn_adv_pra.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/ICE/icedyn_adv_pra.F90

    r12252 r12340  
    4747   !! * Substitutions 
    4848#  include "vectopt_loop_substitute.h90" 
     49#  include "do_loop_substitute.h90" 
    4950   !!---------------------------------------------------------------------- 
    5051   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    102103      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
    103104      DO jl = 1, jpl 
    104          DO jj = 2, jpjm1 
    105             DO ji = fs_2, fs_jpim1 
    106                zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    107                   &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    108                   &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    109                   &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    110                zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    111                   &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    112                   &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    113                   &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    114                zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    115                   &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    116                   &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    117                   &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    118             END DO 
    119          END DO 
     105         DO_2D_00_00 
     106            zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     107               &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     108               &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     109               &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     110            zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     111               &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     112               &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     113               &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     114            zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     115               &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     116               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     117               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     118         END_2D 
    120119      END DO 
    121120      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
    … …  
    252251         ! derive open water from ice concentration 
    253252         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    254          DO jj = 2, jpjm1 
    255             DO ji = fs_2, fs_jpim1 
    256                pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
    257                   &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    258             END DO 
    259          END DO 
     253         DO_2D_00_00 
     254            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
     255               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     256         END_2D 
    260257         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. ) 
    261258         ! 
    … …  
    309306         ! 
    310307         ! Limitation of moments.                                            
    311          DO jj = 2, jpjm1 
    312             DO ji = 1, jpi 
    313                !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    314                psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    315                ! 
    316                zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    317                zs1max  = 1.5 * zslpmax 
    318                zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    319                zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    320                   &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    321                rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    322  
    323                ps0 (ji,jj,jl) = zslpmax   
    324                psx (ji,jj,jl) = zs1new         * rswitch 
    325                psxx(ji,jj,jl) = zs2new         * rswitch 
    326                psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    327                psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    328                psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    329             END DO 
    330          END DO 
     308         DO_2D_00_11 
     309            !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     310            psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
     311            ! 
     312            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     313            zs1max  = 1.5 * zslpmax 
     314            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
     315            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
     316               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
     317            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     318 
     319            ps0 (ji,jj,jl) = zslpmax   
     320            psx (ji,jj,jl) = zs1new         * rswitch 
     321            psxx(ji,jj,jl) = zs2new         * rswitch 
     322            psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
     323            psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
     324            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     325         END_2D 
    331326 
    332327         !  Calculate fluxes and moments between boxes i<-->i+1               
    333          DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0  
    334             DO ji = 1, jpi 
    335                zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    336                zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    337                zalfq        =  zalf * zalf 
    338                zalf1        =  1.0 - zalf 
    339                zalf1q       =  zalf1 * zalf1 
    340                ! 
    341                zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    342                zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    343                zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    344                zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    345                zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    346                zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    347                zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    348  
    349                !  Readjust moments remaining in the box. 
    350                psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    351                ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    352                psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    353                psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    354                psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    355                psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    356                psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    357             END DO 
    358          END DO 
    359  
    360          DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0. 
    361             DO ji = 1, fs_jpim1 
    362                zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    363                zalg  (ji,jj) = zalf 
    364                zalfq         = zalf * zalf 
    365                zalf1         = 1.0 - zalf 
    366                zalg1 (ji,jj) = zalf1 
    367                zalf1q        = zalf1 * zalf1 
    368                zalg1q(ji,jj) = zalf1q 
    369                ! 
    370                zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
    371                zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
    372                   &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
    373                zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
    374                zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
    375                zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
    376                zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
    377                zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    378             END DO 
    379          END DO 
    380  
    381          DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
    382             DO ji = fs_2, fs_jpim1 
    383                zbt  =       zbet(ji-1,jj) 
    384                zbt1 = 1.0 - zbet(ji-1,jj) 
    385                ! 
    386                psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    387                ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    388                psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    389                psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    390                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    391                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    392                psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    393             END DO 
    394          END DO 
     328         DO_2D_00_11 
     329            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     330            zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
     331            zalfq        =  zalf * zalf 
     332            zalf1        =  1.0 - zalf 
     333            zalf1q       =  zalf1 * zalf1 
     334            ! 
     335            zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
     336            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
     337            zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
     338            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
     339            zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     340            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
     341            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
     342 
     343            !  Readjust moments remaining in the box. 
     344            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     345            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     346            psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
     347            psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
     348            psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
     349            psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
     350            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     351         END_2D 
     352 
     353         DO_2D_00_10 
     354            zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     355            zalg  (ji,jj) = zalf 
     356            zalfq         = zalf * zalf 
     357            zalf1         = 1.0 - zalf 
     358            zalg1 (ji,jj) = zalf1 
     359            zalf1q        = zalf1 * zalf1 
     360            zalg1q(ji,jj) = zalf1q 
     361            ! 
     362            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     363            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     364               &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     365            zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     366            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     367            zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     368            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     369            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     370         END_2D 
     371 
     372         DO_2D_00_00 
     373            zbt  =       zbet(ji-1,jj) 
     374            zbt1 = 1.0 - zbet(ji-1,jj) 
     375            ! 
     376            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
     377            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
     378            psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
     379            psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
     380            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
     381            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
     382            psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
     383         END_2D 
    395384 
    396385         !   Put the temporary moments into appropriate neighboring boxes.     
    397          DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0. 
    398             DO ji = fs_2, fs_jpim1 
    399                zbt  =       zbet(ji-1,jj) 
    400                zbt1 = 1.0 - zbet(ji-1,jj) 
    401                psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    402                zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    403                zalf1         = 1.0 - zalf 
    404                ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    405                ! 
    406                ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    407                psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    408                psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    409                   &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
    410                   &            + zbt1 * psxx(ji,jj,jl) 
    411                psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    412                   &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    413                   &            + zbt1 * psxy(ji,jj,jl) 
    414                psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    415                psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    416             END DO 
    417          END DO 
    418  
    419          DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0. 
    420             DO ji = fs_2, fs_jpim1 
    421                zbt  =       zbet(ji,jj) 
    422                zbt1 = 1.0 - zbet(ji,jj) 
    423                psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    424                zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    425                zalf1         = 1.0 - zalf 
    426                ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    427                ! 
    428                ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    429                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    430                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    431                   &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    432                   &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    433                psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    434                   &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    435                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    436                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
    437             END DO 
    438          END DO 
     386         DO_2D_00_00 
     387            zbt  =       zbet(ji-1,jj) 
     388            zbt1 = 1.0 - zbet(ji-1,jj) 
     389            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
     390            zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
     391            zalf1         = 1.0 - zalf 
     392            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
     393            ! 
     394            ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
     395            psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
     396            psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
     397               &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
     398               &            + zbt1 * psxx(ji,jj,jl) 
     399            psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
     400               &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
     401               &            + zbt1 * psxy(ji,jj,jl) 
     402            psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
     403            psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
     404         END_2D 
     405 
     406         DO_2D_00_00 
     407            zbt  =       zbet(ji,jj) 
     408            zbt1 = 1.0 - zbet(ji,jj) 
     409            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     410            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     411            zalf1         = 1.0 - zalf 
     412            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     413            ! 
     414            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
     415            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
     416            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
     417               &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
     418               &                                           + ( zalf1 - zalf ) * ztemp ) ) 
     419            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     420               &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
     421            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
     422            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
     423         END_2D 
    439424 
    440425      END DO 
    … …  
    478463         ! 
    479464         ! Limitation of moments. 
    480          DO jj = 1, jpj 
    481             DO ji = fs_2, fs_jpim1 
    482                !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    483                psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    484                ! 
    485                zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    486                zs1max  = 1.5 * zslpmax 
    487                zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    488                zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    489                   &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
    490                rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    491                ! 
    492                ps0 (ji,jj,jl) = zslpmax   
    493                psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    494                psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    495                psy (ji,jj,jl) = zs1new         * rswitch 
    496                psyy(ji,jj,jl) = zs2new         * rswitch 
    497                psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    498             END DO 
    499          END DO 
     465         DO_2D_11_00 
     466            !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
     467            psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
     468            ! 
     469            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     470            zs1max  = 1.5 * zslpmax 
     471            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
     472            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
     473               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     474            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     475            ! 
     476            ps0 (ji,jj,jl) = zslpmax   
     477            psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
     478            psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
     479            psy (ji,jj,jl) = zs1new         * rswitch 
     480            psyy(ji,jj,jl) = zs2new         * rswitch 
     481            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     482         END_2D 
    500483  
    501484         !  Calculate fluxes and moments between boxes j<-->j+1               
    502          DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
    503             DO ji = fs_2, fs_jpim1 
    504                zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    505                zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
    506                zalfq        =  zalf * zalf 
    507                zalf1        =  1.0 - zalf 
    508                zalf1q       =  zalf1 * zalf1 
    509                ! 
    510                zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    511                zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    512                zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    513                zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    514                zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    515                zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    516                zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
    517                ! 
    518                !  Readjust moments remaining in the box. 
    519                psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    520                ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    521                psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    522                psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    523                psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    524                psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    525                psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    526             END DO 
    527          END DO 
    528          ! 
    529          DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    530             DO ji = fs_2, fs_jpim1 
    531                zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    532                zalg  (ji,jj) = zalf 
    533                zalfq         = zalf * zalf 
    534                zalf1         = 1.0 - zalf 
    535                zalg1 (ji,jj) = zalf1 
    536                zalf1q        = zalf1 * zalf1 
    537                zalg1q(ji,jj) = zalf1q 
    538                ! 
    539                zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl) 
    540                zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) & 
    541                   &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 
    542                zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 
    543                zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq 
    544                zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 
    545                zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl) 
    546                zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl) 
    547             END DO 
    548          END DO 
     485         DO_2D_11_00 
     486            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
     487            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     488            zalfq        =  zalf * zalf 
     489            zalf1        =  1.0 - zalf 
     490            zalf1q       =  zalf1 * zalf1 
     491            ! 
     492            zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
     493            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
     494            zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
     495            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
     496            zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     497            zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
     498            zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
     499            ! 
     500            !  Readjust moments remaining in the box. 
     501            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     502            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     503            psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
     504            psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
     505            psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
     506            psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
     507            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     508         END_2D 
     509         ! 
     510         DO_2D_10_00 
     511            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
     512            zalg  (ji,jj) = zalf 
     513            zalfq         = zalf * zalf 
     514            zalf1         = 1.0 - zalf 
     515            zalg1 (ji,jj) = zalf1 
     516            zalf1q        = zalf1 * zalf1 
     517            zalg1q(ji,jj) = zalf1q 
     518            ! 
     519            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl) 
     520            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) & 
     521               &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 
     522            zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 
     523            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq 
     524            zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 
     525            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl) 
     526            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl) 
     527         END_2D 
    549528 
    550529         !  Readjust moments remaining in the box.  
    551          DO jj = 2, jpjm1 
    552             DO ji = fs_2, fs_jpim1 
    553                zbt  =         zbet(ji,jj-1) 
    554                zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    555                ! 
    556                psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    557                ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    558                psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    559                psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    560                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    561                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    562                psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
    563             END DO 
    564          END DO 
     530         DO_2D_00_00 
     531            zbt  =         zbet(ji,jj-1) 
     532            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     533            ! 
     534            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
     535            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
     536            psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
     537            psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
     538            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
     539            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
     540            psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     541         END_2D 
    565542 
    566543         !   Put the temporary moments into appropriate neighboring boxes.     
    567          DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
    568             DO ji = fs_2, fs_jpim1 
    569                zbt  =       zbet(ji,jj-1) 
    570                zbt1 = 1.0 - zbet(ji,jj-1) 
    571                psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    572                zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    573                zalf1         = 1.0 - zalf 
    574                ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    575                ! 
    576                ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    577                psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    578                   &             + zbt1 * psy(ji,jj,jl)   
    579                psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    580                   &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    581                   &             + zbt1 * psyy(ji,jj,jl) 
    582                psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    583                   &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    584                   &             + zbt1 * psxy(ji,jj,jl) 
    585                psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    586                psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
    587             END DO 
    588          END DO 
    589  
    590          DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0. 
    591             DO ji = fs_2, fs_jpim1 
    592                zbt  =       zbet(ji,jj) 
    593                zbt1 = 1.0 - zbet(ji,jj) 
    594                psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    595                zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    596                zalf1         = 1.0 - zalf 
    597                ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    598                ! 
    599                ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    600                psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    601                psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    602                   &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    603                   &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    604                psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    605                   &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    606                psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    607                psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
    608             END DO 
    609          END DO 
     544         DO_2D_00_00 
     545            zbt  =       zbet(ji,jj-1) 
     546            zbt1 = 1.0 - zbet(ji,jj-1) 
     547            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
     548            zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
     549            zalf1         = 1.0 - zalf 
     550            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
     551            ! 
     552            ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
     553            psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
     554               &             + zbt1 * psy(ji,jj,jl)   
     555            psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
     556               &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     557               &             + zbt1 * psyy(ji,jj,jl) 
     558            psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
     559               &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
     560               &             + zbt1 * psxy(ji,jj,jl) 
     561            psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
     562            psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
     563         END_2D 
     564 
     565         DO_2D_00_00 
     566            zbt  =       zbet(ji,jj) 
     567            zbt1 = 1.0 - zbet(ji,jj) 
     568            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     569            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     570            zalf1         = 1.0 - zalf 
     571            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     572            ! 
     573            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
     574            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
     575            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
     576               &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
     577               &                                            + ( zalf1 - zalf ) * ztemp ) ) 
     578            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     579               &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
     580            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
     581            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
     582         END_2D 
    610583 
    611584      END DO 
    … …  
    646619      DO jl = 1, jpl 
    647620 
    648          DO jj = 1, jpj 
    649             DO ji = 1, jpi 
    650                IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     621         DO_2D_11_11 
     622            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     623               ! 
     624               !                               ! -- check h_ip -- ! 
     625               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     626               IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     627                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     628                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     629                     pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     630                  ENDIF 
     631               ENDIF 
     632               ! 
     633               !                               ! -- check h_i -- ! 
     634               ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     635               zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     636               IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     637                  pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     638               ENDIF 
     639               ! 
     640               !                               ! -- check h_s -- ! 
     641               ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     642               zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     643               IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     644                  zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    651645                  ! 
    652                   !                               ! -- check h_ip -- ! 
    653                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    654                   IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    655                      zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    656                      IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    657                         pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
    658                      ENDIF 
    659                   ENDIF 
     646                  wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
     647                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    660648                  ! 
    661                   !                               ! -- check h_i -- ! 
    662                   ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
    663                   zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
    664                   IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
    665                      pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
    666                   ENDIF 
    667                   ! 
    668                   !                               ! -- check h_s -- ! 
    669                   ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
    670                   zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
    671                   IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
    672                      zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    673                      ! 
    674                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
    675                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    676                      ! 
    677                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
    678                      pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    679                   ENDIF            
    680                   !                   
    681                ENDIF 
    682             END DO 
    683          END DO 
     649                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     650                  pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     651               ENDIF            
     652               !                   
     653            ENDIF 
     654         END_2D 
    684655      END DO  
    685656      ! 
    … …  
    714685      ! -- check snow load -- ! 
    715686      DO jl = 1, jpl 
    716          DO jj = 1, jpj 
    717             DO ji = 1, jpi 
    718                IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    719                   ! 
    720                   zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    721                   ! 
    722                   IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
    723                      ! put snow excess in the ocean 
    724                      zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
    725                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
    726                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    727                      ! correct snow volume and heat content 
    728                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
    729                      pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
    730                   ENDIF 
    731                   ! 
     687         DO_2D_11_11 
     688            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     689               ! 
     690               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     691               ! 
     692               IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
     693                  ! put snow excess in the ocean 
     694                  zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     695                  wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     696                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     697                  ! correct snow volume and heat content 
     698                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     699                  pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
    732700               ENDIF 
    733             END DO 
    734          END DO 
     701               ! 
     702            ENDIF 
     703         END_2D 
    735704      END DO 
    736705      ! 
Note: See TracChangeset for help on using the changeset viewer.