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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icedyn_adv_pra.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icedyn_adv_pra.F90

    r10425 r13463  
    1616   !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero 
    1717   !!---------------------------------------------------------------------- 
     18   USE phycst         ! physical constant 
    1819   USE dom_oce        ! ocean domain 
    1920   USE ice            ! sea-ice variables 
    2021   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call 
     22   USE icevar         ! sea-ice: operations 
    2123   ! 
    2224   USE in_out_manager ! I/O manager 
     
    2527   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    2628   USE lbclnk         ! lateral boundary conditions (or mpp links) 
    27    USE prtctl         ! Print control 
    2829 
    2930   IMPLICIT NONE 
     
    3637   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   ! ice thickness  
    3738   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness 
    38    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! lead fraction 
     39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! ice concentration 
    3940   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity 
    4041   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice 
    4242   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content 
    4343   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content 
     
    4646 
    4747   !! * Substitutions 
    48 #  include "vectopt_loop_substitute.h90" 
     48#  include "do_loop_substitute.h90" 
    4949   !!---------------------------------------------------------------------- 
    5050   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    5454CONTAINS 
    5555 
    56    SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
     56   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    5757      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    5858      !!---------------------------------------------------------------------- 
     
    7070      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    7171      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     72      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     73      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     74      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    7275      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    7376      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     
    8184      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8285      ! 
    83       INTEGER  ::   jk, jl, jt              ! dummy loop indices 
    84       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    85       REAL(wp) ::   zcfl , zusnit           !   -      - 
    86       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea 
    87       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw 
    88       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    89       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp 
    90       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0es 
    91       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei 
     86      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     87      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     88      REAL(wp) ::   zdt                     !   -      - 
     89      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     90      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     91      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
     92      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
     93      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
     95      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
     97      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
    9298      !!---------------------------------------------------------------------- 
    9399      ! 
    94100      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    95101      ! 
    96       ALLOCATE( zarea(jpi,jpj)    , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) ,                       & 
    97          &      z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   & 
    98          &      z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) ) 
    99       ! 
    100       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !         
    101       zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
    102       zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    103       CALL mpp_max( 'icedyn_adv_pra', zcfl ) 
     102      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     103      DO jl = 1, jpl 
     104         DO_2D( 0, 0, 0, 0 ) 
     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 
     118      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 ) 
     120      ! 
     121      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     122      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     123      !              this should not affect too much the stability 
     124      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) ) 
     125      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) ) 
    104126       
    105       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    106       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     127      ! non-blocking global communication send zcflnow and receive zcflprv 
     128      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 
     129 
     130      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2 
     131      ELSE                         ;   icycle = 1 
    107132      ENDIF 
     133      zdt = rDt_ice / REAL(icycle) 
    108134       
    109       zarea(:,:) = e1e2t(:,:) 
    110       !------------------------- 
    111       ! transported fields                                         
    112       !------------------------- 
    113       z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area  
    114       DO jl = 1, jpl 
    115          z0snw(:,:,jl) = pv_s (:,:,  jl) * e1e2t(:,:)     ! Snow volume 
    116          z0ice(:,:,jl) = pv_i (:,:,  jl) * e1e2t(:,:)     ! Ice  volume 
    117          z0ai (:,:,jl) = pa_i (:,:,  jl) * e1e2t(:,:)     ! Ice area 
    118          z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content 
    119          z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content 
    120          DO jk = 1, nlay_s 
    121             z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
     135      ! --- transport --- ! 
     136      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
     137      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
     138 
     139      DO jt = 1, icycle 
     140 
     141         ! record at_i before advection (for open water) 
     142         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     143          
     144         ! --- transported fields --- !                                         
     145         DO jl = 1, jpl 
     146            zarea(:,:,jl) = e1e2t(:,:) 
     147            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume 
     148            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume 
     149            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area 
     150            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content 
     151            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content 
     152            DO jk = 1, nlay_s 
     153               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
     154            END DO 
     155            DO jk = 1, nlay_i 
     156               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     157            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 
     161            ENDIF 
    122162         END DO 
    123          DO jk = 1, nlay_i 
    124             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     163         ! 
     164         !                                                                  !--------------------------------------------! 
     165         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     166            !                                                               !--------------------------------------------! 
     167            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     168            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     169            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     170            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     171            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     172            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     173            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     174            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     175            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     176            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     177            ! 
     178            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     179               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     180                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     181               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     182                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     183            END DO 
     184            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     185               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     186                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     187               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     188                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     189            END DO 
     190            ! 
     191            IF ( ln_pnd_H12 ) THEN 
     192               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     193               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
     194               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     195               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     196            ENDIF 
     197            !                                                               !--------------------------------------------! 
     198         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==! 
     199            !                                                               !--------------------------------------------! 
     200            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     201            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     202            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     203            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     204            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     205            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     206            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     207            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     208            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     209            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     210            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     211               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     212                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     213               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     214                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     215            END DO 
     216            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     217               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     218                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     219               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     220                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     221            END DO 
     222            IF ( ln_pnd_H12 ) THEN 
     223               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     224               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
     225               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     226               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     227            ENDIF 
     228            ! 
     229         ENDIF 
     230 
     231         ! --- Recover the properties from their contents --- ! 
     232         DO jl = 1, jpl 
     233            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     234            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     235            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     236            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     237            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     238            DO jk = 1, nlay_s 
     239               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     240            END DO 
     241            DO jk = 1, nlay_i 
     242               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     243            END DO 
     244            IF ( ln_pnd_H12 ) THEN 
     245               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     246               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     247            ENDIF 
    125248         END DO 
    126          IF ( ln_pnd_H12 ) THEN 
    127             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    128             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
    129          ENDIF 
     249         ! 
     250         ! derive open water from ice concentration 
     251         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     252         DO_2D( 0, 0, 0, 0 ) 
     253            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
     254               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     255         END_2D 
     256         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp ) 
     257         ! 
     258         ! --- Ensure non-negative fields --- ! 
     259         !     Remove negative values (conservation is ensured) 
     260         !     (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 ) 
     262         ! 
     263         ! --- Make sure ice thickness is not too big --- ! 
     264         !     (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 ) 
     266         ! 
     267         ! --- Ensure snow load is not too big --- ! 
     268         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     269         ! 
    130270      END DO 
    131  
    132       !                                                    !--------------------------------------------! 
    133       IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    134          !                                                 !--------------------------------------------! 
    135          DO jt = 1, initad 
    136             CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    137                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    138             CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    139                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    140             DO jl = 1, jpl 
    141                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    142                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    143                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    144                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    145                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    146                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    147                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    148                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    149                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    150                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    151                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    152                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    153                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    154                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    155                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    156                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    157                CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    158                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    159                CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    160                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    161                DO jk = 1, nlay_s                                                               !--- snow heat contents --- 
    162                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    163                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    164                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    165                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    166                END DO 
    167                DO jk = 1, nlay_i                                                               !--- ice heat contents --- 
    168                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    169                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    170                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    171                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    172                END DO 
    173                IF ( ln_pnd_H12 ) THEN 
    174                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
    175                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    176                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &  
    177                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    178                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   -- 
    179                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    180                   CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &  
    181                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    182                ENDIF 
    183             END DO 
    184          END DO 
    185       !                                                    !--------------------------------------------! 
    186       ELSE                                                 !== even ice time step:  adv_y then adv_x  ==! 
    187          !                                                 !--------------------------------------------! 
    188          DO jt = 1, initad 
    189             CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    190                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    191             CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   & 
    192                &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    193             DO jl = 1, jpl 
    194                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    195                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    196                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   & 
    197                   &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    198                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    199                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    200                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    201                   &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    202                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    203                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    204                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    205                   &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    206                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    207                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    208                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   & 
    209                   &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    210                CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    211                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    212                CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    213                   &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    214                DO jk = 1, nlay_s                                                             !--- snow heat contents --- 
    215                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    216                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    217                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   & 
    218                      &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) ) 
    219                END DO 
    220                DO jk = 1, nlay_i                                                             !--- ice heat contents --- 
    221                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    222                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    223                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   &  
    224                      &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    225                END DO 
    226                IF ( ln_pnd_H12 ) THEN 
    227                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
    228                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    229                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
    230                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    231                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   --- 
    232                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    233                   CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
    234                      &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  ) 
    235                ENDIF 
    236             END DO 
    237          END DO 
    238       ENDIF 
    239  
    240       !------------------------------------------- 
    241       ! Recover the properties from their contents 
    242       !------------------------------------------- 
    243       pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1) 
    244       DO jl = 1, jpl 
    245          pv_i (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    246          pv_s (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    247          psv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    248          poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    249          pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    250          DO jk = 1, nlay_s 
    251             pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    252          END DO 
    253          DO jk = 1, nlay_i 
    254             pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    255          END DO 
    256          IF ( ln_pnd_H12 ) THEN 
    257             pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    258             pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    259          ENDIF 
    260       END DO 
    261       ! 
    262       DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei ) 
    263271      ! 
    264272      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file 
     
    267275    
    268276    
    269    SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   & 
     277   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   & 
    270278      &              psx, psxx, psy , psyy, psxy ) 
    271279      !!---------------------------------------------------------------------- 
     
    275283      !!                variable on x axis 
    276284      !!---------------------------------------------------------------------- 
    277       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    278       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    279       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
    280       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
    281       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
    282       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
    283       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     285      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step 
     286      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     287      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
     288      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     289      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     290      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     291      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    284292      !!  
    285       INTEGER  ::   ji, jj                               ! dummy loop indices 
    286       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars 
     293      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     294      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    287295      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    288296      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     
    291299      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    292300      !----------------------------------------------------------------------- 
    293  
    294       ! Limitation of moments.                                            
    295  
    296       zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection. 
    297  
    298       DO jj = 1, jpj 
    299          DO ji = 1, jpi 
    300             zslpmax = MAX( 0._wp, ps0(ji,jj) ) 
     301      ! 
     302      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     303      ! 
     304      DO jl = 1, jcat   ! loop on categories 
     305         ! 
     306         ! Limitation of moments.                                            
     307         DO_2D( 0, 0, 1, 1 ) 
     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) ) 
    301312            zs1max  = 1.5 * zslpmax 
    302             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) ) 
     313            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    303314            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    304                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
     315               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    305316            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    306317 
    307             ps0 (ji,jj) = zslpmax   
    308             psx (ji,jj) = zs1new      * rswitch 
    309             psxx(ji,jj) = zs2new      * rswitch 
    310             psy (ji,jj) = psy (ji,jj) * rswitch 
    311             psyy(ji,jj) = psyy(ji,jj) * rswitch 
    312             psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 
    313          END DO 
    314       END DO 
    315  
    316       !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    317       psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    318  
    319       !  Calculate fluxes and moments between boxes i<-->i+1               
    320       DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0  
    321          DO ji = 1, jpi 
     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( 0, 0, 1, 1 ) 
    322328            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    323             zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj) 
     329            zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    324330            zalfq        =  zalf * zalf 
    325331            zalf1        =  1.0 - zalf 
    326332            zalf1q       =  zalf1 * zalf1 
    327333            ! 
    328             zfm (ji,jj)  =  zalf  *   psm (ji,jj) 
    329             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) ) 
    330             zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) ) 
    331             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq 
    332             zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) ) 
    333             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj) 
    334             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj) 
     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) 
    335341 
    336342            !  Readjust moments remaining in the box. 
    337             psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    338             ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    339             psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) ) 
    340             psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj) 
    341             psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj) 
    342             psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj) 
    343             psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    344          END DO 
    345       END DO 
    346  
    347       DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0. 
    348          DO ji = 1, fs_jpim1 
    349             zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj)  
     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( 0, 0, 1, 0 ) 
     353            zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    350354            zalg  (ji,jj) = zalf 
    351355            zalfq         = zalf * zalf 
     
    355359            zalg1q(ji,jj) = zalf1q 
    356360            ! 
    357             zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj) 
    358             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) ) 
    359             zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) ) 
    360             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq 
    361             zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) ) 
    362             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj) 
    363             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj) 
    364          END DO 
    365       END DO 
    366  
    367       DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
    368          DO ji = fs_2, fs_jpim1 
     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( 0, 0, 0, 0 ) 
    369372            zbt  =       zbet(ji-1,jj) 
    370373            zbt1 = 1.0 - zbet(ji-1,jj) 
    371374            ! 
    372             psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) ) 
    373             ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) ) 
    374             psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) ) 
    375             psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj) 
    376             psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) ) 
    377             psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) ) 
    378             psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj) 
    379          END DO 
    380       END DO 
    381  
    382       !   Put the temporary moments into appropriate neighboring boxes.     
    383       DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0. 
    384          DO ji = fs_2, fs_jpim1 
     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( 0, 0, 0, 0 ) 
    385386            zbt  =       zbet(ji-1,jj) 
    386387            zbt1 = 1.0 - zbet(ji-1,jj) 
    387             psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj) 
    388             zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj) 
    389             zalf1       = 1.0 - zalf 
    390             ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj) 
    391             ! 
    392             ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj) 
    393             psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj) 
    394             psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               & 
    395                &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   & 
    396                &                                                + zbt1 * psxx(ji,jj) 
    397             psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             & 
    398                &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   & 
    399                &                                                + zbt1 * psxy(ji,jj) 
    400             psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj) 
    401             psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj) 
    402          END DO 
    403       END DO 
    404  
    405       DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0. 
    406          DO ji = fs_2, fs_jpim1 
     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( 0, 0, 0, 0 ) 
    407406            zbt  =       zbet(ji,jj) 
    408407            zbt1 = 1.0 - zbet(ji,jj) 
    409             psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    410             zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    411             zalf1       = 1.0 - zalf 
    412             ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    413             ! 
    414             ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    415             psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) 
    416             psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  & 
    417                &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      & 
    418                &                                      + ( zalf1 - zalf ) * ztemp ) ) 
    419             psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  & 
    420                &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  ) 
    421             psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) ) 
    422             psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) ) 
    423          END DO 
     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 
    424424      END DO 
    425425 
    426426      !-- Lateral boundary conditions 
    427       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1., ps0 , 'T',  1.   & 
    428          &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
    429          &              , psxx, 'T',  1., psyy, 'T',  1.   & 
    430          &              , psxy, 'T',  1. ) 
    431  
    432       IF(ln_ctl) THEN 
    433          CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    434          CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    435          CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    436          CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 
    437       ENDIF 
     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 ) 
    438430      ! 
    439431   END SUBROUTINE adv_x 
    440432 
    441433 
    442    SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   & 
     434   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   & 
    443435      &              psx, psxx, psy , psyy, psxy ) 
    444436      !!--------------------------------------------------------------------- 
     
    448440      !!                variable on y axis 
    449441      !!--------------------------------------------------------------------- 
    450       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    451       REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
    452       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
    453       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area 
    454       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected 
    455       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments  
    456       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    457       !! 
    458       INTEGER  ::   ji, jj                               ! dummy loop indices 
    459       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars 
     442      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step 
     443      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     444      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
     445      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     446      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     447      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     448      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
     449      !! 
     450      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     451      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    460452      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    461453      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     
    464456      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    465457      !--------------------------------------------------------------------- 
    466  
    467       ! Limitation of moments. 
    468  
    469       zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 
    470  
    471       DO jj = 1, jpj 
    472          DO ji = 1, jpi 
    473             zslpmax = MAX( 0._wp, ps0(ji,jj) ) 
     458      ! 
     459      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     460      !       
     461      DO jl = 1, jcat   ! loop on categories 
     462         ! 
     463         ! Limitation of moments. 
     464         DO_2D( 1, 1, 0, 0 ) 
     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) ) 
    474469            zs1max  = 1.5 * zslpmax 
    475             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
     470            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    476471            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    477                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
     472               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
    478473            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    479474            ! 
    480             ps0 (ji,jj) = zslpmax   
    481             psx (ji,jj) = psx (ji,jj) * rswitch 
    482             psxx(ji,jj) = psxx(ji,jj) * rswitch 
    483             psy (ji,jj) = zs1new * rswitch 
    484             psyy(ji,jj) = zs2new * rswitch 
    485             psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 
    486          END DO 
    487       END DO 
    488  
    489       !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    490       psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  ) 
    491  
    492       !  Calculate fluxes and moments between boxes j<-->j+1               
    493       DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
    494          DO ji = 1, jpi 
     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( 1, 1, 0, 0 ) 
    495485            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    496             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
     486            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
    497487            zalfq        =  zalf * zalf 
    498488            zalf1        =  1.0 - zalf 
    499489            zalf1q       =  zalf1 * zalf1 
    500490            ! 
    501             zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    502             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj) ) )  
    503             zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
    504             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
    505             zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
    506             zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    507             zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
     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) 
    508498            ! 
    509499            !  Readjust moments remaining in the box. 
    510             psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    511             ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    512             psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
    513             psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
    514             psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
    515             psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
    516             psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    517          END DO 
    518       END DO 
    519       ! 
    520       DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    521          DO ji = 1, jpi 
    522             zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
     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) 
     507         END_2D 
     508         ! 
     509         DO_2D( 1, 0, 0, 0 ) 
     510            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    523511            zalg  (ji,jj) = zalf 
    524512            zalfq         = zalf * zalf 
     
    528516            zalg1q(ji,jj) = zalf1q 
    529517            ! 
    530             zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1) 
    531             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
    532             zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
    533             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq 
    534             zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
    535             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1) 
    536             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1) 
    537          END DO 
    538       END DO 
    539  
    540       !  Readjust moments remaining in the box.  
    541       DO jj = 2, jpj 
    542          DO ji = 1, jpi 
     518            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl) 
     519            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) & 
     520               &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 
     521            zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 
     522            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq 
     523            zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 
     524            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl) 
     525            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl) 
     526         END_2D 
     527 
     528         !  Readjust moments remaining in the box.  
     529         DO_2D( 0, 0, 0, 0 ) 
    543530            zbt  =         zbet(ji,jj-1) 
    544531            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    545532            ! 
    546             psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
    547             ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
    548             psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
    549             psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 
    550             psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
    551             psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
    552             psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
    553          END DO 
     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) 
     540         END_2D 
     541 
     542         !   Put the temporary moments into appropriate neighboring boxes.     
     543         DO_2D( 0, 0, 0, 0 ) 
     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( 0, 0, 0, 0 ) 
     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 
    554583      END DO 
    555584 
    556       !   Put the temporary moments into appropriate neighboring boxes.     
    557       DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
    558          DO ji = 1, jpi 
    559             zbt  =         zbet(ji,jj-1) 
    560             zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    561             psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
    562             zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
    563             zalf1       = 1.0 - zalf 
    564             ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
    565             ! 
    566             ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj) 
    567             psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
    568                &                                               + zbt1 * psy(ji,jj)   
    569             psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
    570                &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
    571                &                                               + zbt1 * psyy(ji,jj) 
    572             psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               & 
    573                &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   & 
    574                &                                                + zbt1 * psxy(ji,jj) 
    575             psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
    576             psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
    577          END DO 
     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 ) 
     589      ! 
     590   END SUBROUTINE adv_y 
     591 
     592 
     593   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     594      !!------------------------------------------------------------------- 
     595      !!                  ***  ROUTINE Hbig  *** 
     596      !! 
     597      !! ** Purpose : Thickness correction in case advection scheme creates 
     598      !!              abnormally tick ice or snow 
     599      !! 
     600      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     601      !!                 (before advection) and reduce it by adapting ice concentration 
     602      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     603      !!                 (before advection) and reduce it by sending the excess in the ocean 
     604      !! 
     605      !! ** input   : Max thickness of the surrounding 9-points 
     606      !!------------------------------------------------------------------- 
     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 
     610      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 
     614      !!------------------------------------------------------------------- 
     615      ! 
     616      z1_dt = 1._wp / pdt 
     617      ! 
     618      DO jl = 1, jpl 
     619 
     620         DO_2D( 1, 1, 1, 1 ) 
     621            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     622               ! 
     623               !                               ! -- check h_ip -- ! 
     624               ! 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 
     626                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     627                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     628                     pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     629                  ENDIF 
     630               ENDIF 
     631               ! 
     632               !                               ! -- check h_i -- ! 
     633               ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     634               zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     635               IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     636                  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) 
     637               ENDIF 
     638               ! 
     639               !                               ! -- check h_s -- ! 
     640               ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     641               zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     642               IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     643                  zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     644                  ! 
     645                  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 
     646                  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 
     647                  ! 
     648                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     649                  pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     650               ENDIF            
     651               !                   
     652            ENDIF 
     653         END_2D 
     654      END DO  
     655      ! 
     656   END SUBROUTINE Hbig 
     657 
     658 
     659   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     660      !!------------------------------------------------------------------- 
     661      !!                  ***  ROUTINE Hsnow  *** 
     662      !! 
     663      !! ** Purpose : 1- Check snow load after advection 
     664      !!              2- Correct pond concentration to avoid a_ip > a_i 
     665      !! 
     666      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface 
     667      !!              then put the snow excess in the ocean 
     668      !! 
     669      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards 
     670      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 
     671      !!              make the snow very thick (if concentration decreases drastically) 
     672      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 
     673      !!------------------------------------------------------------------- 
     674      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step 
     675      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip 
     676      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     677      ! 
     678      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     679      REAL(wp) ::   z1_dt, zvs_excess, zfra 
     680      !!------------------------------------------------------------------- 
     681      ! 
     682      z1_dt = 1._wp / pdt 
     683      ! 
     684      ! -- check snow load -- ! 
     685      DO jl = 1, jpl 
     686         DO_2D( 1, 1, 1, 1 ) 
     687            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     688               ! 
     689               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos ) 
     690               ! 
     691               IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
     692                  ! put snow excess in the ocean 
     693                  zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     694                  wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     695                  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 
     696                  ! correct snow volume and heat content 
     697                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     698                  pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
     699               ENDIF 
     700               ! 
     701            ENDIF 
     702         END_2D 
    578703      END DO 
    579  
    580       DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
    581          DO ji = 1, jpi 
    582             zbt  =         zbet(ji,jj) 
    583             zbt1 = ( 1.0 - zbet(ji,jj) ) 
    584             psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    585             zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    586             zalf1       = 1.0 - zalf 
    587             ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj) 
    588             ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) ) 
    589             psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
    590             psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
    591                &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          & 
    592                &                                         + ( zalf1 - zalf ) * ztemp )                                ) 
    593             psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       & 
    594                &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  ) 
    595             psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
    596             psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
    597          END DO 
    598       END DO 
    599  
    600       !-- Lateral boundary conditions 
    601       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1.,  ps0 , 'T',  1.   & 
    602          &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
    603          &              , psxx, 'T',  1.,  psyy, 'T',  1.   & 
    604          &              , psxy, 'T',  1. ) 
    605  
    606       IF(ln_ctl) THEN 
    607          CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    608          CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    609          CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    610          CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 
    611       ENDIF 
    612       ! 
    613    END SUBROUTINE adv_y 
     704      ! 
     705      !-- correct pond concentration to avoid a_ip > a_i -- ! 
     706      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
     707      ! 
     708   END SUBROUTINE Hsnow 
    614709 
    615710 
     
    624719      ! 
    625720      !                             !* allocate prather fields 
    626       ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   & 
    627          &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
     721      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
    628722         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   & 
    629723         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   & 
     
    652746      !!                   ***  ROUTINE adv_pra_rst  *** 
    653747      !!                      
    654       !! ** Purpose :   Read or write RHG file in restart file 
     748      !! ** Purpose :   Read or write file in restart file 
    655749      !! 
    656750      !! ** Method  :   use of IOM library 
     
    671765         !                                   !==========================! 
    672766         ! 
    673          IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0 
     767         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0 
    674768         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0 
    675769         ENDIF 
     
    678772            ! 
    679773            !                                                        ! ice thickness 
    680             CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  ) 
    681             CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  ) 
    682             CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice ) 
    683             CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice ) 
    684             CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice ) 
     774            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice  ) 
     775            CALL iom_get( numrir, jpdom_auto, 'syice' , syice  ) 
     776            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 
     777            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 
     778            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 
    685779            !                                                        ! snow thickness 
    686             CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   ) 
    687             CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   ) 
    688             CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  ) 
    689             CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  ) 
    690             CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  ) 
    691             !                                                        ! lead fraction 
    692             CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    ) 
    693             CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    ) 
    694             CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   ) 
    695             CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   ) 
    696             CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   ) 
     780            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn   ) 
     781            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn   ) 
     782            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  ) 
     783            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  ) 
     784            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  ) 
     785            !                                                        ! ice concentration 
     786            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    ) 
     787            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    ) 
     788            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   ) 
     789            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   ) 
     790            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   ) 
    697791            !                                                        ! ice salinity 
    698             CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal  ) 
    699             CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal  ) 
    700             CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal ) 
    701             CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal ) 
    702             CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal ) 
     792            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal  ) 
     793            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal  ) 
     794            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 
     795            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 
     796            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 
    703797            !                                                        ! ice age 
    704             CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage  ) 
    705             CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage  ) 
    706             CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage ) 
    707             CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 
    708             CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 
    709             !                                                        ! open water in sea ice 
    710             CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw  ) 
    711             CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw  ) 
    712             CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw ) 
    713             CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw ) 
    714             CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw ) 
     798            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage  ) 
     799            CALL iom_get( numrir, jpdom_auto, 'syage' , syage  ) 
     800            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 
     801            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 
     802            CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage ) 
    715803            !                                                        ! snow layers heat content 
    716804            DO jk = 1, nlay_s 
    717805               WRITE(zchar1,'(I2.2)') jk 
    718                znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
    719                znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
    720                znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
    721                znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
    722                znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
     806               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
     807               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
     808               znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
     809               znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
     810               znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
    723811            END DO 
    724812            !                                                        ! ice layers heat content 
    725813            DO jk = 1, nlay_i 
    726814               WRITE(zchar1,'(I2.2)') jk 
    727                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
    728                znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
    729                znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
    730                znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
    731                znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_autoglo, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
     815               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
     816               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
     817               znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
     818               znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
     819               znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
    732820            END DO 
    733821            ! 
    734822            IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction 
    735                CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
    736                CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  ) 
    737                CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap ) 
    738                CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap ) 
    739                CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap ) 
     823               CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  ) 
     824               CALL iom_get( numrir, jpdom_auto, 'syap' , syap  ) 
     825               CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
     826               CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
     827               CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
    740828               !                                                     ! melt pond volume 
    741                CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  ) 
    742                CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  ) 
    743                CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp ) 
    744                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
    745                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     829               CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  ) 
     830               CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  ) 
     831               CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
     832               CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
     833               CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
    746834            ENDIF 
    747835            ! 
     
    752840            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness 
    753841            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness 
    754             sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction 
     842            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration 
    755843            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity 
    756844            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age 
    757             sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice 
    758845            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    759846            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
     
    786873         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  ) 
    787874         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  ) 
    788          !                                                           ! lead fraction 
     875         !                                                           ! ice concentration 
    789876         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    ) 
    790877         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    ) 
     
    804891         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 
    805892         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 
    806          !                                                           ! open water in sea ice 
    807          CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw  ) 
    808          CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw  ) 
    809          CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw ) 
    810          CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw ) 
    811          CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw ) 
    812893         !                                                           ! snow layers heat content 
    813894         DO jk = 1, nlay_s 
Note: See TracChangeset for help on using the changeset viewer.