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 13899 for NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_adv_pra.F90 – NEMO

Ignore:
Timestamp:
2020-11-27T17:26:33+01:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: update branch to trunk and add ICB test case

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/tickets_icb_1900

    • Property svn:externals

      •  

        old new  
        22^/utils/build/makenemo@HEAD   makenemo 
        33^/utils/build/mk@HEAD         mk 
        4 ^/utils/tools/@HEAD           tools 
         4^/utils/tools@HEAD            tools 
        55^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        … …  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@12931        sette 
         10^/utils/CI/sette@13559        sette 

  • NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_adv_pra.F90

    r13226 r13899  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume 
    4647 
    4748   !! * Substitutions 
    … …  
    5556 
    5657   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
    … …  
    8182      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    8283      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness 
    8385      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8486      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8587      ! 
    86       INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     88      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8789      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    88       REAL(wp) ::   zdt                     !   -      - 
     90      REAL(wp) ::   zdt, z1_dt              !   -      - 
    8991      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9092      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    9193      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
    92       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     95      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    9397      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9498      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    95       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl 
    96100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    97101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     102      !! diagnostics 
     103      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    98104      !!---------------------------------------------------------------------- 
    99105      ! 
    100106      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    101107      ! 
    102       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
    103       DO jl = 1, jpl 
    104          DO_2D_00_00 
    105             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    106                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    107                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    108                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    109             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    110                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    111                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    112                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    113             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    114                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    115                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    116                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    117          END_2D 
     108      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     109      ! thickness and salinity 
     110      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     112      END WHERE 
     113      CALL icemax3D( ph_i , zhi_max ) 
     114      CALL icemax3D( ph_s , zhs_max ) 
     115      CALL icemax3D( ph_ip, zhip_max) 
     116      CALL icemax3D( zs_i , zsi_max ) 
     117      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
     118      ! 
     119      ! enthalpies 
     120      DO jk = 1, nlay_i 
     121         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     122         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     123         END WHERE 
    118124      END DO 
    119       CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 
     125      DO jk = 1, nlay_s 
     126         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     128         END WHERE 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
     132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
     134      ! 
    120135      ! 
    121136      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
    … …  
    132147      ENDIF 
    133148      zdt = rDt_ice / REAL(icycle) 
     149      z1_dt = 1._wp / zdt 
    134150       
    135151      ! --- transport --- ! 
    … …  
    138154 
    139155      DO jt = 1, icycle 
     156 
     157         ! diagnostics 
     158         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     159         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     160         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     161            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    140162 
    141163         ! record at_i before advection (for open water) 
    … …  
    156178               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    157179            END DO 
    158             IF ( ln_pnd_H12 ) THEN 
    159                z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    160                z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     180            IF ( ln_pnd_LEV ) THEN 
     181               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
     182               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
     183               IF ( ln_pnd_lids ) THEN 
     184                  z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume 
     185               ENDIF 
    161186            ENDIF 
    162187         END DO 
    … …  
    189214            END DO 
    190215            ! 
    191             IF ( ln_pnd_H12 ) THEN 
     216            IF ( ln_pnd_LEV ) THEN 
    192217               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    193218               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    194219               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    195220               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     221               IF ( ln_pnd_lids ) THEN 
     222                  CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     223                  CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     224               ENDIF 
    196225            ENDIF 
    197226            !                                                               !--------------------------------------------! 
    … …  
    220249                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    221250            END DO 
    222             IF ( ln_pnd_H12 ) THEN 
     251            IF ( ln_pnd_LEV ) THEN 
    223252               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    224253               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    225254               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    226255               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     256               IF ( ln_pnd_lids ) THEN 
     257                  CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     258                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     259               ENDIF 
    227260            ENDIF 
    228261            ! 
     262         ENDIF 
     263          
     264         ! --- Lateral boundary conditions --- ! 
     265         !     caution: for gradients (sx and sy) the sign changes 
     266         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume 
     267            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  & 
     268            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume 
     269            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  ) 
     270         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity 
     271            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  & 
     272            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration 
     273            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  ) 
     274         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age 
     275            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  ) 
     276         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy 
     277            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  )  
     278         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
     279            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
     280         IF ( ln_pnd_LEV ) THEN 
     281            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
     282               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
     283               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume 
     284               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )  
     285            IF ( ln_pnd_lids ) THEN 
     286               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume 
     287                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )  
     288            ENDIF 
    229289         ENDIF 
    230290 
    … …  
    242302               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    243303            END DO 
    244             IF ( ln_pnd_H12 ) THEN 
     304            IF ( ln_pnd_LEV ) THEN 
    245305               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    246306               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     307               IF ( ln_pnd_lids ) THEN 
     308                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     309               ENDIF 
    247310            ENDIF 
    248311         END DO 
    … …  
    250313         ! derive open water from ice concentration 
    251314         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    252          DO_2D_00_00 
     315         DO_2D( 0, 0, 0, 0 ) 
    253316            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
    254317               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    … …  
    256319         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp ) 
    257320         ! 
     321         ! --- diagnostics --- ! 
     322         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     323            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     324         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     325            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     326         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     327            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     328            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
     329         ! 
    258330         ! --- Ensure non-negative fields --- ! 
    259331         !     Remove negative values (conservation is ensured) 
    260332         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    261          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     333         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    262334         ! 
    263335         ! --- Make sure ice thickness is not too big --- ! 
    264336         !     (because ice thickness can be too large where ice concentration is very small) 
    265          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     337         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     338            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    266339         ! 
    267340         ! --- Ensure snow load is not too big --- ! 
    … …  
    292365      !!  
    293366      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     367      INTEGER  ::   jj0                                  ! dummy loop indices 
    294368      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    295369      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    296370      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     371      REAL(wp) ::   zpsm, zps0 
     372      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    297373      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    298374      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    299375      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    300376      !----------------------------------------------------------------------- 
     377      ! in order to avoid lbc_lnk (communications): 
     378      !    jj loop must be 1:jpj   if adv_x is called first 
     379      !                and 2:jpj-1 if adv_x is called second 
     380      jj0 = NINT(pcrh) 
    301381      ! 
    302382      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    305385         ! 
    306386         ! Limitation of moments.                                            
    307          DO_2D_00_11 
    308             !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    309             psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    310             ! 
    311             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    312             zs1max  = 1.5 * zslpmax 
    313             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    314             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    315                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    316             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    317  
    318             ps0 (ji,jj,jl) = zslpmax   
    319             psx (ji,jj,jl) = zs1new         * rswitch 
    320             psxx(ji,jj,jl) = zs2new         * rswitch 
    321             psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    322             psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    323             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    324          END_2D 
    325  
    326          !  Calculate fluxes and moments between boxes i<-->i+1               
    327          DO_2D_00_11 
    328             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    329             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    330             zalfq        =  zalf * zalf 
    331             zalf1        =  1.0 - zalf 
    332             zalf1q       =  zalf1 * zalf1 
    333             ! 
    334             zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    335             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    336             zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    337             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    338             zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    339             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    340             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    341  
    342             !  Readjust moments remaining in the box. 
    343             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    344             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    345             psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    346             psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    347             psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    348             psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    349             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    350          END_2D 
    351  
    352          DO_2D_00_10 
    353             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    354             zalg  (ji,jj) = zalf 
    355             zalfq         = zalf * zalf 
    356             zalf1         = 1.0 - zalf 
    357             zalg1 (ji,jj) = zalf1 
    358             zalf1q        = zalf1 * zalf1 
    359             zalg1q(ji,jj) = zalf1q 
    360             ! 
    361             zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
    362             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
    363                &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
    364             zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
    365             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
    366             zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
    367             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
    368             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    369          END_2D 
    370  
    371          DO_2D_00_00 
    372             zbt  =       zbet(ji-1,jj) 
    373             zbt1 = 1.0 - zbet(ji-1,jj) 
    374             ! 
    375             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    376             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    377             psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    378             psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    379             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    380             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    381             psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    382          END_2D 
    383  
    384          !   Put the temporary moments into appropriate neighboring boxes.     
    385          DO_2D_00_00 
    386             zbt  =       zbet(ji-1,jj) 
    387             zbt1 = 1.0 - zbet(ji-1,jj) 
    388             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    389             zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    390             zalf1         = 1.0 - zalf 
    391             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    392             ! 
    393             ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    394             psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    395             psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    396                &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
    397                &            + zbt1 * psxx(ji,jj,jl) 
    398             psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    399                &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    400                &            + zbt1 * psxy(ji,jj,jl) 
    401             psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    402             psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    403          END_2D 
    404  
    405          DO_2D_00_00 
    406             zbt  =       zbet(ji,jj) 
    407             zbt1 = 1.0 - zbet(ji,jj) 
    408             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    409             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    410             zalf1         = 1.0 - zalf 
    411             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    412             ! 
    413             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    414             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    415             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    416                &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    417                &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    418             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    419                &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    420             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    421             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
    422          END_2D 
    423  
     387         DO jj = Njs0 - jj0, Nje0 + jj0 
     388             
     389            DO ji = Nis0 - 1, Nie0 + 1 
     390 
     391               zpsm  = psm (ji,jj,jl) ! optimization 
     392               zps0  = ps0 (ji,jj,jl) 
     393               zpsx  = psx (ji,jj,jl) 
     394               zpsxx = psxx(ji,jj,jl) 
     395               zpsy  = psy (ji,jj,jl) 
     396               zpsyy = psyy(ji,jj,jl) 
     397               zpsxy = psxy(ji,jj,jl) 
     398 
     399               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     400               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     401               ! 
     402               zslpmax = MAX( 0._wp, zps0 ) 
     403               zs1max  = 1.5 * zslpmax 
     404               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     405               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
     406               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     407 
     408               zps0  = zslpmax   
     409               zpsx  = zs1new  * rswitch 
     410               zpsxx = zs2new  * rswitch 
     411               zpsy  = zpsy    * rswitch 
     412               zpsyy = zpsyy   * rswitch 
     413               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     414 
     415               !  Calculate fluxes and moments between boxes i<-->i+1               
     416               !                                !  Flux from i to i+1 WHEN u GT 0  
     417               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     418               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
     419               zalfq        =  zalf * zalf 
     420               zalf1        =  1.0 - zalf 
     421               zalf1q       =  zalf1 * zalf1 
     422               ! 
     423               zfm (ji,jj)  =  zalf  *   zpsm  
     424               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     425               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     426               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     427               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     428               zfxy(ji,jj)  =  zalfq *   zpsxy 
     429               zfyy(ji,jj)  =  zalf  *   zpsyy 
     430 
     431               !                                !  Readjust moments remaining in the box. 
     432               zpsm  =  zpsm  - zfm(ji,jj) 
     433               zps0  =  zps0  - zf0(ji,jj) 
     434               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     435               zpsxx =  zalf1  * zalf1q * zpsxx 
     436               zpsy  =  zpsy  - zfy (ji,jj) 
     437               zpsyy =  zpsyy - zfyy(ji,jj) 
     438               zpsxy =  zalf1q * zpsxy 
     439               ! 
     440               psm (ji,jj,jl) = zpsm ! optimization 
     441               ps0 (ji,jj,jl) = zps0  
     442               psx (ji,jj,jl) = zpsx  
     443               psxx(ji,jj,jl) = zpsxx 
     444               psy (ji,jj,jl) = zpsy  
     445               psyy(ji,jj,jl) = zpsyy 
     446               psxy(ji,jj,jl) = zpsxy 
     447               ! 
     448            END DO 
     449             
     450            DO ji = Nis0 - 1, Nie0 
     451               !                                !  Flux from i+1 to i when u LT 0. 
     452               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     453               zalg  (ji,jj) = zalf 
     454               zalfq         = zalf * zalf 
     455               zalf1         = 1.0 - zalf 
     456               zalg1 (ji,jj) = zalf1 
     457               zalf1q        = zalf1 * zalf1 
     458               zalg1q(ji,jj) = zalf1q 
     459               ! 
     460               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     461               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     462                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     463               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     464               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     465               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     466               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     467               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     468            END DO 
     469 
     470            DO ji = Nis0, Nie0 
     471               ! 
     472               zpsm  = psm (ji,jj,jl) ! optimization 
     473               zps0  = ps0 (ji,jj,jl) 
     474               zpsx  = psx (ji,jj,jl) 
     475               zpsxx = psxx(ji,jj,jl) 
     476               zpsy  = psy (ji,jj,jl) 
     477               zpsyy = psyy(ji,jj,jl) 
     478               zpsxy = psxy(ji,jj,jl) 
     479               !                                !  Readjust moments remaining in the box. 
     480               zbt  =       zbet(ji-1,jj) 
     481               zbt1 = 1.0 - zbet(ji-1,jj) 
     482               ! 
     483               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     484               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     485               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     486               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     487               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) ) 
     488               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     489               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     490 
     491               !   Put the temporary moments into appropriate neighboring boxes.     
     492               !                                !   Flux from i to i+1 IF u GT 0. 
     493               zbt   =       zbet(ji-1,jj) 
     494               zbt1  = 1.0 - zbet(ji-1,jj) 
     495               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     496               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     497               zalf1 = 1.0 - zalf 
     498               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     499               ! 
     500               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     501               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     502               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     503                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     504                  &            + zbt1 * zpsxx 
     505               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     506                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     507                  &            + zbt1 * zpsxy 
     508               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     509               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     510 
     511               !                                !  Flux from i+1 to i IF u LT 0. 
     512               zbt   =       zbet(ji,jj) 
     513               zbt1  = 1.0 - zbet(ji,jj) 
     514               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     515               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     516               zalf1 = 1.0 - zalf 
     517               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     518               ! 
     519               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     520               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     521               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     522                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     523                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     524               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     525                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     526               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     527               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     528               ! 
     529               psm (ji,jj,jl) = zpsm  ! optimization 
     530               ps0 (ji,jj,jl) = zps0  
     531               psx (ji,jj,jl) = zpsx  
     532               psxx(ji,jj,jl) = zpsxx 
     533               psy (ji,jj,jl) = zpsy  
     534               psyy(ji,jj,jl) = zpsyy 
     535               psxy(ji,jj,jl) = zpsxy 
     536            END DO 
     537            ! 
     538         END DO 
     539         ! 
    424540      END DO 
    425  
    426       !-- Lateral boundary conditions 
    427       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    428          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    429          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    430       ! 
     541      !       
    431542   END SUBROUTINE adv_x 
    432543 
    … …  
    449560      !! 
    450561      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     562      INTEGER  ::   ji0                                  ! dummy loop indices 
    451563      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    452564      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    453565      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     566      REAL(wp) ::   zpsm, zps0 
     567      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    454568      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    455569      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    456570      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    457571      !--------------------------------------------------------------------- 
     572      ! in order to avoid lbc_lnk (communications): 
     573      !    ji loop must be 1:jpi   if adv_y is called first 
     574      !                and 2:jpi-1 if adv_y is called second 
     575      ji0 = NINT(pcrh) 
    458576      ! 
    459577      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    462580         ! 
    463581         ! Limitation of moments. 
    464          DO_2D_11_00 
    465             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    466             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    467             ! 
    468             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     582         DO_2D( 1, 1, ji0, ji0 ) 
     583            ! 
     584            zpsm  = psm (ji,jj,jl) ! optimization 
     585            zps0  = ps0 (ji,jj,jl) 
     586            zpsx  = psx (ji,jj,jl) 
     587            zpsxx = psxx(ji,jj,jl) 
     588            zpsy  = psy (ji,jj,jl) 
     589            zpsyy = psyy(ji,jj,jl) 
     590            zpsxy = psxy(ji,jj,jl) 
     591            ! 
     592            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     593            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     594            ! 
     595            zslpmax = MAX( 0._wp, zps0 ) 
    469596            zs1max  = 1.5 * zslpmax 
    470             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    471             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    472                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     597            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     598            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    473599            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    474600            ! 
    475             ps0 (ji,jj,jl) = zslpmax   
    476             psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    477             psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    478             psy (ji,jj,jl) = zs1new         * rswitch 
    479             psyy(ji,jj,jl) = zs2new         * rswitch 
    480             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    481          END_2D 
    482   
    483          !  Calculate fluxes and moments between boxes j<-->j+1               
    484          DO_2D_11_00 
     601            zps0  = zslpmax   
     602            zpsx  = zpsx  * rswitch 
     603            zpsxx = zpsxx * rswitch 
     604            zpsy  = zs1new         * rswitch 
     605            zpsyy = zs2new         * rswitch 
     606            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     607 
     608            !  Calculate fluxes and moments between boxes j<-->j+1               
     609            !                                !  Flux from j to j+1 WHEN v GT 0    
    485610            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    486             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     611            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    487612            zalfq        =  zalf * zalf 
    488613            zalf1        =  1.0 - zalf 
    489614            zalf1q       =  zalf1 * zalf1 
    490615            ! 
    491             zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    492             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    493             zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    494             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    495             zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    496             zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    497             zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
    498             ! 
    499             !  Readjust moments remaining in the box. 
    500             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    501             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    502             psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    503             psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    504             psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    505             psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    506             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     616            zfm (ji,jj)  =  zalf  * zpsm 
     617            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     618            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     619            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     620            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     621            zfxy(ji,jj)  =  zalfq * zpsxy 
     622            zfxx(ji,jj)  =  zalf  * zpsxx 
     623            ! 
     624            !                                !  Readjust moments remaining in the box. 
     625            zpsm   =  zpsm  - zfm(ji,jj) 
     626            zps0   =  zps0  - zf0(ji,jj) 
     627            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     628            zpsyy  =  zalf1 * zalf1q * zpsyy 
     629            zpsx   =  zpsx  - zfx(ji,jj) 
     630            zpsxx  =  zpsxx - zfxx(ji,jj) 
     631            zpsxy  =  zalf1q * zpsxy 
     632            ! 
     633            psm (ji,jj,jl) = zpsm ! optimization 
     634            ps0 (ji,jj,jl) = zps0  
     635            psx (ji,jj,jl) = zpsx  
     636            psxx(ji,jj,jl) = zpsxx 
     637            psy (ji,jj,jl) = zpsy  
     638            psyy(ji,jj,jl) = zpsyy 
     639            psxy(ji,jj,jl) = zpsxy 
    507640         END_2D 
    508641         ! 
    509          DO_2D_10_00 
     642         DO_2D( 1, 0, ji0, ji0 ) 
     643            !                                !  Flux from j+1 to j when v LT 0. 
    510644            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    511645            zalg  (ji,jj) = zalf 
    … …  
    526660         END_2D 
    527661 
    528          !  Readjust moments remaining in the box.  
    529          DO_2D_00_00 
     662         DO_2D( 0, 0, ji0, ji0 ) 
     663            !                                !  Readjust moments remaining in the box. 
    530664            zbt  =         zbet(ji,jj-1) 
    531665            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    532666            ! 
    533             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    534             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    535             psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    536             psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    537             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    538             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    539             psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     667            zpsm  = psm (ji,jj,jl) ! optimization 
     668            zps0  = ps0 (ji,jj,jl) 
     669            zpsx  = psx (ji,jj,jl) 
     670            zpsxx = psxx(ji,jj,jl) 
     671            zpsy  = psy (ji,jj,jl) 
     672            zpsyy = psyy(ji,jj,jl) 
     673            zpsxy = psxy(ji,jj,jl) 
     674            ! 
     675            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     676            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     677            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     678            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     679            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     680            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     681            zpsxy = zalg1q(ji,jj-1) * zpsxy 
     682 
     683            !   Put the temporary moments into appropriate neighboring boxes.     
     684            !                                !   Flux from j to j+1 IF v GT 0. 
     685            zbt   =       zbet(ji,jj-1) 
     686            zbt1  = 1.0 - zbet(ji,jj-1) 
     687            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     688            zalf  = zbt * zfm(ji,jj-1) / zpsm  
     689            zalf1 = 1.0 - zalf 
     690            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     691            ! 
     692            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     693            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     694               &             + zbt1 * zpsy   
     695            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     696               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     697               &             + zbt1 * zpsyy 
     698            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     699               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     700               &             + zbt1 * zpsxy 
     701            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     702            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
     703 
     704            !                                !  Flux from j+1 to j IF v LT 0. 
     705            zbt   =       zbet(ji,jj) 
     706            zbt1  = 1.0 - zbet(ji,jj) 
     707            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     708            zalf  = zbt1 * zfm(ji,jj) / zpsm 
     709            zalf1 = 1.0 - zalf 
     710            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     711            ! 
     712            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     713            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     714            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     715               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     716               &                         + ( zalf1 - zalf ) * ztemp ) ) 
     717            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     718               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     719            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     720            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     721            ! 
     722            psm (ji,jj,jl) = zpsm ! optimization 
     723            ps0 (ji,jj,jl) = zps0  
     724            psx (ji,jj,jl) = zpsx  
     725            psxx(ji,jj,jl) = zpsxx 
     726            psy (ji,jj,jl) = zpsy  
     727            psyy(ji,jj,jl) = zpsyy 
     728            psxy(ji,jj,jl) = zpsxy 
    540729         END_2D 
    541  
    542          !   Put the temporary moments into appropriate neighboring boxes.     
    543          DO_2D_00_00 
    544             zbt  =       zbet(ji,jj-1) 
    545             zbt1 = 1.0 - zbet(ji,jj-1) 
    546             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    547             zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    548             zalf1         = 1.0 - zalf 
    549             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    550             ! 
    551             ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    552             psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    553                &             + zbt1 * psy(ji,jj,jl)   
    554             psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    555                &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    556                &             + zbt1 * psyy(ji,jj,jl) 
    557             psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    558                &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    559                &             + zbt1 * psxy(ji,jj,jl) 
    560             psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    561             psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
    562          END_2D 
    563  
    564          DO_2D_00_00 
    565             zbt  =       zbet(ji,jj) 
    566             zbt1 = 1.0 - zbet(ji,jj) 
    567             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    568             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    569             zalf1         = 1.0 - zalf 
    570             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    571             ! 
    572             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    573             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    574             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    575                &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    576                &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    577             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    578                &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    579             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    580             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
    581          END_2D 
    582  
     730         ! 
    583731      END DO 
    584  
    585       !-- Lateral boundary conditions 
    586       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    587          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    588          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    589732      ! 
    590733   END SUBROUTINE adv_y 
    591734 
    592735 
    593    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     736   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     737      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    594738      !!------------------------------------------------------------------- 
    595739      !!                  ***  ROUTINE Hbig  *** 
    … …  
    605749      !! ** input   : Max thickness of the surrounding 9-points 
    606750      !!------------------------------------------------------------------- 
    607       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    608       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    609       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     751      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     752      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     753      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     754      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     755      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    610756      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    611       ! 
    612       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    613       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     757      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     758      ! 
     759      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     760      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    614761      !!------------------------------------------------------------------- 
    615762      ! 
    … …  
    617764      ! 
    618765      DO jl = 1, jpl 
    619  
    620          DO_2D_11_11 
     766         DO_2D( 1, 1, 1, 1 ) 
    621767            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    622768               ! 
    623769               !                               ! -- check h_ip -- ! 
    624770               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    625                IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     771               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    626772                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    627773                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    … …  
    650796               ENDIF            
    651797               !                   
     798               !                               ! -- check s_i -- ! 
     799               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     800               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     801               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     802                  zfra = psi_max(ji,jj,jl) / zsi 
     803                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     804                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     805               ENDIF 
     806               ! 
    652807            ENDIF 
    653808         END_2D 
    654809      END DO  
     810      ! 
     811      !                                           ! -- check e_i/v_i -- ! 
     812      DO jl = 1, jpl 
     813         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     814            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     815               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     816               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     817               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     818                  zfra = pei_max(ji,jj,jk,jl) / zei 
     819                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     820                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     821               ENDIF 
     822            ENDIF 
     823         END_3D 
     824      END DO 
     825      !                                           ! -- check e_s/v_s -- ! 
     826      DO jl = 1, jpl 
     827         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     828            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     829               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     830               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     831               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     832                  zfra = pes_max(ji,jj,jk,jl) / zes 
     833                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     834                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     835               ENDIF 
     836            ENDIF 
     837         END_3D 
     838      END DO 
    655839      ! 
    656840   END SUBROUTINE Hbig 
    … …  
    684868      ! -- check snow load -- ! 
    685869      DO jl = 1, jpl 
    686          DO_2D_11_11 
     870         DO_2D( 1, 1, 1, 1 ) 
    687871            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    688872               ! 
    … …  
    724908         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   & 
    725909         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   & 
    726          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    727          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     910         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
     911         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     912         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   & 
    728913         ! 
    729914         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
    … …  
    772957            ! 
    773958            !                                                        ! ice thickness 
    774             CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice ) 
    775             CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice ) 
    776             CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice ) 
    777             CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice ) 
    778             CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice ) 
     959            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 
     960            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 
     961            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 
     962            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 
     963            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 
    779964            !                                                        ! snow thickness 
    780             CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn  ) 
    781             CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn  ) 
    782             CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  ) 
    783             CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  ) 
    784             CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  ) 
     965            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp ) 
     966            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp ) 
     967            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  ) 
     968            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  ) 
     969            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  ) 
    785970            !                                                        ! ice concentration 
    786             CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    ) 
    787             CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    ) 
    788             CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   ) 
    789             CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   ) 
    790             CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   ) 
     971            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp ) 
     972            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp ) 
     973            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   ) 
     974            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   ) 
     975            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   ) 
    791976            !                                                        ! ice salinity 
    792             CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal ) 
    793             CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal ) 
    794             CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal ) 
    795             CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal ) 
    796             CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal ) 
     977            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 
     978            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 
     979            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 
     980            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 
     981            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 
    797982            !                                                        ! ice age 
    798             CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage ) 
    799             CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage ) 
    800             CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage ) 
    801             CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 
    802             CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 
     983            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 
     984            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 
     985            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 
     986            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 
     987            CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage ) 
    803988            !                                                        ! snow layers heat content 
    804989            DO jk = 1, nlay_s 
    805990               WRITE(zchar1,'(I2.2)') jk 
    806                znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
    807                znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
    808                znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
    809                znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
    810                znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
     991               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
     992               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
     993               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
     994               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
     995               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
    811996            END DO 
    812997            !                                                        ! ice layers heat content 
    813998            DO jk = 1, nlay_i 
    814999               WRITE(zchar1,'(I2.2)') jk 
    815                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
    816                znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
    817                znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
    818                znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
    819                znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
    820             END DO 
    821             ! 
    822             IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction 
    823                CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
    824                CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  ) 
    825                CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap ) 
    826                CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap ) 
    827                CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap ) 
    828                !                                                     ! melt pond volume 
    829                CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  ) 
    830                CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  ) 
    831                CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp ) 
    832                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
    833                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     1000               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
     1001               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
     1002               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
     1003               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
     1004               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
     1005            END DO 
     1006            ! 
     1007            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     1008               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     1009                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 
     1010                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 
     1011                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
     1012                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
     1013                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
     1014                  !                                                     ! melt pond volume 
     1015                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 
     1016                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 
     1017                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
     1018                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
     1019                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
     1020               ELSE 
     1021                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction 
     1022                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume 
     1023               ENDIF 
     1024                  ! 
     1025               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
     1026                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     1027                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 
     1028                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 
     1029                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 
     1030                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 
     1031                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl ) 
     1032                  ELSE 
     1033                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume 
     1034                  ENDIF 
     1035               ENDIF 
    8341036            ENDIF 
    8351037            ! 
    … …  
    8451047            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    8461048            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    847             IF( ln_pnd_H12 ) THEN 
    848                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    849                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     1049            IF( ln_pnd_LEV ) THEN 
     1050               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
     1051               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
     1052               IF ( ln_pnd_lids ) THEN 
     1053                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume 
     1054               ENDIF 
    8501055            ENDIF 
    8511056         ENDIF 
    … …  
    9101115         END DO 
    9111116         ! 
    912          IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction 
     1117         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
    9131118            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    9141119            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 
    … …  
    9221127            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    9231128            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 
     1129            ! 
     1130            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume 
     1131               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  ) 
     1132               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  ) 
     1133               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl ) 
     1134               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl ) 
     1135               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl ) 
     1136            ENDIF 
    9241137         ENDIF 
    9251138         ! 
    … …  
    9271140      ! 
    9281141   END SUBROUTINE adv_pra_rst 
     1142 
     1143   SUBROUTINE icemax3D( pice , pmax ) 
     1144      !!--------------------------------------------------------------------- 
     1145      !!                   ***  ROUTINE icemax3D ***                      
     1146      !! ** Purpose :  compute the max of the 9 points around 
     1147      !!---------------------------------------------------------------------- 
     1148      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1149      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1150      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1151      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1152      !!---------------------------------------------------------------------- 
     1153      DO jl = 1, jpl 
     1154         DO jj = Njs0-1, Nje0+1     
     1155            DO ji = Nis0, Nie0 
     1156               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1157            END DO 
     1158         END DO 
     1159         DO jj = Njs0, Nje0     
     1160            DO ji = Nis0, Nie0 
     1161               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1162            END DO 
     1163         END DO 
     1164      END DO 
     1165   END SUBROUTINE icemax3D 
     1166 
     1167   SUBROUTINE icemax4D( pice , pmax ) 
     1168      !!--------------------------------------------------------------------- 
     1169      !!                   ***  ROUTINE icemax4D ***                      
     1170      !! ** Purpose :  compute the max of the 9 points around 
     1171      !!---------------------------------------------------------------------- 
     1172      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1173      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1174      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1175      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1176      !!---------------------------------------------------------------------- 
     1177      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1178      DO jl = 1, jpl 
     1179         DO jk = 1, jlay 
     1180            DO jj = Njs0-1, Nje0+1     
     1181               DO ji = Nis0, Nie0 
     1182                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1183               END DO 
     1184            END DO 
     1185            DO jj = Njs0, Nje0     
     1186               DO ji = Nis0, Nie0 
     1187                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1188               END DO 
     1189            END DO 
     1190         END DO 
     1191      END DO 
     1192   END SUBROUTINE icemax4D 
    9291193 
    9301194#else 
Note: See TracChangeset for help on using the changeset viewer.