Ignore:
Timestamp:
2019-12-11T09:27:27+01:00 (10 months ago)
Author:
davestorkey
Message:

2019/dev_ASINTER-01-05_merged: Update to r12072 of trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_pra.F90

    r10425 r12165  
    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 
     
    8181      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8282      ! 
    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 
     83      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     84      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
     85      REAL(wp) ::   zdt                     !   -      - 
     86      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     87      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     88      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
     89      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
     90      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
     91      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     92      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
     93      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
    9294      !!---------------------------------------------------------------------- 
    9395      ! 
    9496      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    9597      ! 
    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 ) 
     98      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     99      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     100      !              this should not affect too much the stability 
     101      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
     102      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    104103       
    105       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    106       ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     104      ! non-blocking global communication send zcflnow and receive zcflprv 
     105      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 
     106 
     107      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2 
     108      ELSE                         ;   icycle = 1 
    107109      ENDIF 
     110      zdt = rdt_ice / REAL(icycle) 
    108111       
    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 
    122          END DO 
    123          DO jk = 1, nlay_i 
    124             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    125          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 
     112      ! --- transport --- ! 
     113      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
     114      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
     115 
     116      DO jt = 1, icycle 
     117 
     118         ! record at_i before advection (for open water) 
     119         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     120          
     121         ! --- transported fields --- !                                         
     122         DO jl = 1, jpl 
     123            zarea(:,:,jl) = e1e2t(:,:) 
     124            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume 
     125            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume 
     126            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area 
     127            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content 
     128            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content 
     129            DO jk = 1, nlay_s 
     130               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
     131            END DO 
     132            DO jk = 1, nlay_i 
     133               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     134            END DO 
     135            IF ( ln_pnd_H12 ) THEN 
     136               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
     137               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     138            ENDIF 
     139         END DO 
     140         ! 
     141         !                                                                  !--------------------------------------------! 
     142         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     143            !                                                               !--------------------------------------------! 
     144            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     145            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     146            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     147            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     148            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     149            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     150            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     151            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     152            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     153            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     154            ! 
     155            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     156               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     157                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     158               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     159                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     160            END DO 
     161            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     162               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     163                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     164               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     165                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     166            END DO 
     167            ! 
     168            IF ( ln_pnd_H12 ) THEN 
     169               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     170               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
     171               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     172               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     173            ENDIF 
     174            !                                                               !--------------------------------------------! 
     175         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==! 
     176            !                                                               !--------------------------------------------! 
     177            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     178            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     179            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     180            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     181            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     182            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     183            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     184            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     185            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     186            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     187            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     188               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     189                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     190               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     191                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     192            END DO 
     193            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     194               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     195                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     196               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     197                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     198            END DO 
     199            IF ( ln_pnd_H12 ) THEN 
     200               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     201               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
     202               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     203               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     204            ENDIF 
     205            ! 
    129206         ENDIF 
     207 
     208         ! --- Recover the properties from their contents --- ! 
     209         DO jl = 1, jpl 
     210            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     211            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     212            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     213            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     214            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     215            DO jk = 1, nlay_s 
     216               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     217            END DO 
     218            DO jk = 1, nlay_i 
     219               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     220            END DO 
     221            IF ( ln_pnd_H12 ) THEN 
     222               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     223               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     224            ENDIF 
     225         END DO 
     226         ! 
     227         ! derive open water from ice concentration 
     228         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     229         DO jj = 2, jpjm1 
     230            DO ji = fs_2, fs_jpim1 
     231               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
     232                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     233            END DO 
     234         END DO 
     235         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. ) 
     236         ! 
     237         ! --- Ensure non-negative fields --- ! 
     238         !     Remove negative values (conservation is ensured) 
     239         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     240         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 ) 
     241         ! 
     242         ! --- Ensure snow load is not too big --- ! 
     243         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     244         ! 
    130245      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 ) 
    263246      ! 
    264247      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file 
     
    267250    
    268251    
    269    SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   & 
     252   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   & 
    270253      &              psx, psxx, psy , psyy, psxy ) 
    271254      !!---------------------------------------------------------------------- 
     
    275258      !!                variable on x axis 
    276259      !!---------------------------------------------------------------------- 
    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 
     260      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step 
     261      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     262      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s] 
     263      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     264      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     265      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     266      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    284267      !!  
    285       INTEGER  ::   ji, jj                               ! dummy loop indices 
    286       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars 
     268      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     269      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    287270      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    288271      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     
    291274      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    292275      !----------------------------------------------------------------------- 
    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             zs1max  = 1.5 * zslpmax 
    302             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) ) 
    303             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    304                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
    305             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    306  
    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 
     276      ! 
     277      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     278      ! 
     279      DO jl = 1, jcat   ! loop on categories 
     280         ! 
     281         ! Limitation of moments.                                            
     282         DO jj = 2, jpjm1 
     283            DO ji = 1, jpi 
     284               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     285               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
     286               ! 
     287               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     288               zs1max  = 1.5 * zslpmax 
     289               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
     290               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
     291                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
     292               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     293 
     294               ps0 (ji,jj,jl) = zslpmax   
     295               psx (ji,jj,jl) = zs1new         * rswitch 
     296               psxx(ji,jj,jl) = zs2new         * rswitch 
     297               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
     298               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
     299               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     300            END DO 
     301         END DO 
     302 
     303         !  Calculate fluxes and moments between boxes i<-->i+1               
     304         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0  
     305            DO ji = 1, jpi 
     306               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     307               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
     308               zalfq        =  zalf * zalf 
     309               zalf1        =  1.0 - zalf 
     310               zalf1q       =  zalf1 * zalf1 
     311               ! 
     312               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
     313               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
     314               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
     315               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
     316               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     317               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
     318               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
     319 
     320               !  Readjust moments remaining in the box. 
     321               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     322               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     323               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
     324               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
     325               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
     326               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
     327               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     328            END DO 
     329         END DO 
     330 
     331         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0. 
     332            DO ji = 1, fs_jpim1 
     333               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     334               zalg  (ji,jj) = zalf 
     335               zalfq         = zalf * zalf 
     336               zalf1         = 1.0 - zalf 
     337               zalg1 (ji,jj) = zalf1 
     338               zalf1q        = zalf1 * zalf1 
     339               zalg1q(ji,jj) = zalf1q 
     340               ! 
     341               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     342               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     343                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     344               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     345               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     346               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     347               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     348               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     349            END DO 
     350         END DO 
     351 
     352         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
     353            DO ji = fs_2, fs_jpim1 
     354               zbt  =       zbet(ji-1,jj) 
     355               zbt1 = 1.0 - zbet(ji-1,jj) 
     356               ! 
     357               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
     358               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
     359               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
     360               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
     361               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
     362               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
     363               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
     364            END DO 
     365         END DO 
     366 
     367         !   Put the temporary moments into appropriate neighboring boxes.     
     368         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0. 
     369            DO ji = fs_2, fs_jpim1 
     370               zbt  =       zbet(ji-1,jj) 
     371               zbt1 = 1.0 - zbet(ji-1,jj) 
     372               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
     373               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
     374               zalf1         = 1.0 - zalf 
     375               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
     376               ! 
     377               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
     378               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
     379               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
     380                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
     381                  &            + zbt1 * psxx(ji,jj,jl) 
     382               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
     383                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
     384                  &            + zbt1 * psxy(ji,jj,jl) 
     385               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
     386               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
     387            END DO 
     388         END DO 
     389 
     390         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0. 
     391            DO ji = fs_2, fs_jpim1 
     392               zbt  =       zbet(ji,jj) 
     393               zbt1 = 1.0 - zbet(ji,jj) 
     394               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     395               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     396               zalf1         = 1.0 - zalf 
     397               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     398               ! 
     399               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
     400               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
     401               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
     402                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
     403                  &                                           + ( zalf1 - zalf ) * ztemp ) ) 
     404               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     405                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
     406               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
     407               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
     408            END DO 
     409         END DO 
     410 
    314411      END DO 
    315412 
    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 
    322             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) 
    324             zalfq        =  zalf * zalf 
    325             zalf1        =  1.0 - zalf 
    326             zalf1q       =  zalf1 * zalf1 
    327             ! 
    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) 
    335  
    336             !  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)  
    350             zalg  (ji,jj) = zalf 
    351             zalfq         = zalf * zalf 
    352             zalf1         = 1.0 - zalf 
    353             zalg1 (ji,jj) = zalf1 
    354             zalf1q        = zalf1 * zalf1 
    355             zalg1q(ji,jj) = zalf1q 
    356             ! 
    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 
    369             zbt  =       zbet(ji-1,jj) 
    370             zbt1 = 1.0 - zbet(ji-1,jj) 
    371             ! 
    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 
    385             zbt  =       zbet(ji-1,jj) 
    386             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 
    407             zbt  =       zbet(ji,jj) 
    408             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 
    424       END DO 
    425  
    426413      !-- 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 
     414      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   & 
     415         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     416         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. ) 
    438417      ! 
    439418   END SUBROUTINE adv_x 
    440419 
    441420 
    442    SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   & 
     421   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   & 
    443422      &              psx, psxx, psy , psyy, psxy ) 
    444423      !!--------------------------------------------------------------------- 
     
    448427      !!                variable on y axis 
    449428      !!--------------------------------------------------------------------- 
    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 
     429      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step 
     430      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0) 
     431      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s] 
     432      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area 
     433      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected 
     434      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments  
     435      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments 
    457436      !! 
    458       INTEGER  ::   ji, jj                               ! dummy loop indices 
    459       REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars 
     437      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     438      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    460439      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    461440      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     
    464443      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    465444      !--------------------------------------------------------------------- 
    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) ) 
    474             zs1max  = 1.5 * zslpmax 
    475             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
    476             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    477                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    478             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    479             ! 
    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 
     445      ! 
     446      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
     447      !       
     448      DO jl = 1, jcat   ! loop on categories 
     449         ! 
     450         ! Limitation of moments. 
     451         DO jj = 1, jpj 
     452            DO ji = fs_2, fs_jpim1 
     453               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
     454               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
     455               ! 
     456               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     457               zs1max  = 1.5 * zslpmax 
     458               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
     459               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
     460                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     461               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     462               ! 
     463               ps0 (ji,jj,jl) = zslpmax   
     464               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
     465               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
     466               psy (ji,jj,jl) = zs1new         * rswitch 
     467               psyy(ji,jj,jl) = zs2new         * rswitch 
     468               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
     469            END DO 
     470         END DO 
     471  
     472         !  Calculate fluxes and moments between boxes j<-->j+1               
     473         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0    
     474            DO ji = fs_2, fs_jpim1 
     475               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
     476               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     477               zalfq        =  zalf * zalf 
     478               zalf1        =  1.0 - zalf 
     479               zalf1q       =  zalf1 * zalf1 
     480               ! 
     481               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
     482               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
     483               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
     484               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
     485               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
     486               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
     487               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
     488               ! 
     489               !  Readjust moments remaining in the box. 
     490               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
     491               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
     492               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
     493               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
     494               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
     495               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
     496               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     497            END DO 
     498         END DO 
     499         ! 
     500         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
     501            DO ji = fs_2, fs_jpim1 
     502               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
     503               zalg  (ji,jj) = zalf 
     504               zalfq         = zalf * zalf 
     505               zalf1         = 1.0 - zalf 
     506               zalg1 (ji,jj) = zalf1 
     507               zalf1q        = zalf1 * zalf1 
     508               zalg1q(ji,jj) = zalf1q 
     509               ! 
     510               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl) 
     511               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) & 
     512                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 
     513               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 
     514               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq 
     515               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 
     516               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl) 
     517               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl) 
     518            END DO 
     519         END DO 
     520 
     521         !  Readjust moments remaining in the box.  
     522         DO jj = 2, jpjm1 
     523            DO ji = fs_2, fs_jpim1 
     524               zbt  =         zbet(ji,jj-1) 
     525               zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     526               ! 
     527               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
     528               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
     529               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
     530               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
     531               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
     532               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
     533               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     534            END DO 
     535         END DO 
     536 
     537         !   Put the temporary moments into appropriate neighboring boxes.     
     538         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
     539            DO ji = fs_2, fs_jpim1 
     540               zbt  =       zbet(ji,jj-1) 
     541               zbt1 = 1.0 - zbet(ji,jj-1) 
     542               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
     543               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
     544               zalf1         = 1.0 - zalf 
     545               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
     546               ! 
     547               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
     548               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
     549                  &             + zbt1 * psy(ji,jj,jl)   
     550               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
     551                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     552                  &             + zbt1 * psyy(ji,jj,jl) 
     553               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
     554                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
     555                  &             + zbt1 * psxy(ji,jj,jl) 
     556               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
     557               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
     558            END DO 
     559         END DO 
     560 
     561         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0. 
     562            DO ji = fs_2, fs_jpim1 
     563               zbt  =       zbet(ji,jj) 
     564               zbt1 = 1.0 - zbet(ji,jj) 
     565               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
     566               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
     567               zalf1         = 1.0 - zalf 
     568               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
     569               ! 
     570               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
     571               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
     572               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
     573                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
     574                  &                                            + ( zalf1 - zalf ) * ztemp ) ) 
     575               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
     576                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
     577               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
     578               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
     579            END DO 
     580         END DO 
     581 
    487582      END DO 
    488583 
    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 
    495             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) 
    497             zalfq        =  zalf * zalf 
    498             zalf1        =  1.0 - zalf 
    499             zalf1q       =  zalf1 * zalf1 
    500             ! 
    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) 
    508             ! 
    509             !  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) 
     584      !-- Lateral boundary conditions 
     585      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   & 
     586         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     587         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. ) 
     588      ! 
     589   END SUBROUTINE adv_y 
     590 
     591 
     592   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     593      !!------------------------------------------------------------------- 
     594      !!                  ***  ROUTINE Hsnow  *** 
     595      !! 
     596      !! ** Purpose : 1- Check snow load after advection 
     597      !!              2- Correct pond concentration to avoid a_ip > a_i 
     598      !! 
     599      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface 
     600      !!              then put the snow excess in the ocean 
     601      !! 
     602      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards 
     603      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 
     604      !!              make the snow very thick (if concentration decreases drastically) 
     605      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 
     606      !!------------------------------------------------------------------- 
     607      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step 
     608      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip 
     609      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     610      ! 
     611      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     612      REAL(wp) ::   z1_dt, zvs_excess, zfra 
     613      !!------------------------------------------------------------------- 
     614      ! 
     615      z1_dt = 1._wp / pdt 
     616      ! 
     617      ! -- check snow load -- ! 
     618      DO jl = 1, jpl 
     619         DO jj = 1, jpj 
     620            DO ji = 1, jpi 
     621               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     622                  ! 
     623                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     624                  ! 
     625                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
     626                     ! put snow excess in the ocean 
     627                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     628                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     629                     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 
     630                     ! correct snow volume and heat content 
     631                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     632                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
     633                  ENDIF 
     634                  ! 
     635               ENDIF 
     636            END DO 
    517637         END DO 
    518638      END DO 
    519639      ! 
    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)  
    523             zalg  (ji,jj) = zalf 
    524             zalfq         = zalf * zalf 
    525             zalf1         = 1.0 - zalf 
    526             zalg1 (ji,jj) = zalf1 
    527             zalf1q        = zalf1 * zalf1 
    528             zalg1q(ji,jj) = zalf1q 
    529             ! 
    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 
    543             zbt  =         zbet(ji,jj-1) 
    544             zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    545             ! 
    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 
    554       END DO 
    555  
    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 
    578       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 
     640      !-- correct pond concentration to avoid a_ip > a_i -- ! 
     641      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
     642      ! 
     643   END SUBROUTINE Hsnow 
    614644 
    615645 
     
    624654      ! 
    625655      !                             !* 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) ,   & 
     656      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
    628657         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   & 
    629658         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   & 
     
    652681      !!                   ***  ROUTINE adv_pra_rst  *** 
    653682      !!                      
    654       !! ** Purpose :   Read or write RHG file in restart file 
     683      !! ** Purpose :   Read or write file in restart file 
    655684      !! 
    656685      !! ** Method  :   use of IOM library 
     
    671700         !                                   !==========================! 
    672701         ! 
    673          IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0 
     702         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0 
    674703         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0 
    675704         ENDIF 
     
    689718            CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  ) 
    690719            CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  ) 
    691             !                                                        ! lead fraction 
     720            !                                                        ! ice concentration 
    692721            CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    ) 
    693722            CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    ) 
     
    707736            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 
    708737            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 ) 
    715738            !                                                        ! snow layers heat content 
    716739            DO jk = 1, nlay_s 
     
    752775            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness 
    753776            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 
     777            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration 
    755778            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity 
    756779            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 
    758780            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    759781            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
     
    786808         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  ) 
    787809         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  ) 
    788          !                                                           ! lead fraction 
     810         !                                                           ! ice concentration 
    789811         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    ) 
    790812         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    ) 
     
    804826         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 
    805827         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 ) 
    812828         !                                                           ! snow layers heat content 
    813829         DO jk = 1, nlay_s 
Note: See TracChangeset for help on using the changeset viewer.