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 13662 for NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_adv_pra.F90 – NEMO

Ignore:
Timestamp:
2020-10-22T20:49:56+02:00 (4 years ago)
Author:
clem
Message:

update to almost r4.0.4

Location:
NEMO/branches/2019/dev_r11842_SI3-10_EAP
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11842_SI3-10_EAP

    • Property svn:externals

      •  

        old new  
        1 ^/utils/build/arch@HEAD       arch 
        2 ^/utils/build/makenemo@HEAD   makenemo 
        3 ^/utils/build/mk@HEAD         mk 
        4 ^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
        6 ^/vendors/FCM@HEAD            ext/FCM 
        7 ^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         1^/utils/build/arch@12130      arch 
         2^/utils/build/makenemo@12191  makenemo 
         3^/utils/build/mk@11662        mk 
         4^/utils/tools_r4.0-HEAD@12672 tools 
         5^/vendors/AGRIF/dev@10586     ext/AGRIF 
         6^/vendors/FCM@10134           ext/FCM 
         7^/vendors/IOIPSL@9655         ext/IOIPSL 
         8 
         9# SETTE mapping (inactive) 
         10#^/utils/CI/sette@12135        sette 

  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_adv_pra.F90

    r11812 r13662  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume 
    4647 
    4748   !! * Substitutions 
    … …  
    5455CONTAINS 
    5556 
    56    SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     57   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
    … …  
    7071      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    7172      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     73      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     74      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     75      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    7276      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    7377      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
    … …  
    7882      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    7983      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness 
    8085      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8186      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8287      ! 
    83       INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     88      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8489      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    85       REAL(wp) ::   zdt                     !   -      - 
     90      REAL(wp) ::   zdt, z1_dt              !   -      - 
    8691      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    8792      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    8893      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     95      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    8997      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9098      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    91       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl 
    92100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    93101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     102      !! diagnostics 
     103      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    94104      !!---------------------------------------------------------------------- 
    95105      ! 
    96106      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
     107      ! 
     108      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     109      ! thickness and salinity 
     110      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     112      END WHERE 
     113      CALL icemax3D( ph_i , zhi_max ) 
     114      CALL icemax3D( ph_s , zhs_max ) 
     115      CALL icemax3D( ph_ip, zhip_max) 
     116      CALL icemax3D( zs_i , zsi_max ) 
     117      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) 
     118 
     119      ! enthalpies 
     120      DO jk = 1, nlay_i 
     121         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     122         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     123         END WHERE 
     124      END DO 
     125      DO jk = 1, nlay_s 
     126         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     128         END WHERE 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
     132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     134      ! 
    97135      ! 
    98136      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
    … …  
    109147      ENDIF 
    110148      zdt = rdt_ice / REAL(icycle) 
     149      z1_dt = 1._wp / zdt 
    111150       
    112151      ! --- transport --- ! 
    … …  
    115154 
    116155      DO jt = 1, icycle 
     156 
     157         ! diagnostics 
     158         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     159         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     160         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     161            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    117162 
    118163         ! record at_i before advection (for open water) 
    … …  
    133178               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    134179            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 
     180            IF ( ln_pnd_LEV ) THEN 
     181               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
     182               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
     183               IF ( ln_pnd_lids ) THEN 
     184                  z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume 
     185               ENDIF 
    138186            ENDIF 
    139187         END DO 
    … …  
    166214            END DO 
    167215            ! 
    168             IF ( ln_pnd_H12 ) THEN 
     216            IF ( ln_pnd_LEV ) THEN 
    169217               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    170218               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    171219               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    172220               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     221               IF ( ln_pnd_lids ) THEN 
     222                  CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     223                  CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     224               ENDIF 
    173225            ENDIF 
    174226            !                                                               !--------------------------------------------! 
    … …  
    197249                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    198250            END DO 
    199             IF ( ln_pnd_H12 ) THEN 
     251            IF ( ln_pnd_LEV ) THEN 
    200252               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    201253               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    202254               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    203255               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     256               IF ( ln_pnd_lids ) THEN 
     257                  CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     258                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     259               ENDIF 
    204260            ENDIF 
    205261            ! 
    206262         ENDIF 
    207  
     263         
     264         ! --- Lateral boundary conditions --- ! 
     265         !     caution: for gradients (sx and sy) the sign changes 
     266         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume 
     267            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  & 
     268            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume 
     269            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  ) 
     270         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity 
     271            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  & 
     272            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration 
     273            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  ) 
     274         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age 
     275            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  ) 
     276         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy 
     277            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  )  
     278         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
     279            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
     280         IF ( ln_pnd_LEV ) THEN 
     281            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
     282               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
     283               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume 
     284               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )  
     285            IF ( ln_pnd_lids ) THEN 
     286               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume 
     287                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )  
     288            ENDIF 
     289         ENDIF 
     290          
    208291         ! --- Recover the properties from their contents --- ! 
    209292         DO jl = 1, jpl 
    … …  
    219302               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    220303            END DO 
    221             IF ( ln_pnd_H12 ) THEN 
     304            IF ( ln_pnd_LEV ) THEN 
    222305               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    223306               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     307               IF ( ln_pnd_lids ) THEN 
     308                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     309               ENDIF 
    224310            ENDIF 
    225311         END DO 
    … …  
    235321         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. ) 
    236322         ! 
     323         ! --- diagnostics --- ! 
     324         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     325            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
     326         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
     327            &                                        - zdiag_adv_salt(:,:) ) * z1_dt 
     328         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     329            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 
     330            &                                        - zdiag_adv_heat(:,:) ) * z1_dt 
     331         ! 
    237332         ! --- Ensure non-negative fields --- ! 
    238333         !     Remove negative values (conservation is ensured) 
    239334         !     (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 ) 
     335         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
     336         ! 
     337         ! --- Make sure ice thickness is not too big --- ! 
     338         !     (because ice thickness can be too large where ice concentration is very small) 
     339         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     340            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    241341         ! 
    242342         ! --- Ensure snow load is not too big --- ! 
    … …  
    267367      !!  
    268368      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     369      INTEGER  ::   jjmin, jjmax                         ! dummy loop indices 
    269370      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    270371      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    271372      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     373      REAL(wp) ::   zpsm, zps0 
     374      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    272375      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    273376      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    274377      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    275378      !----------------------------------------------------------------------- 
     379      ! in order to avoid lbc_lnk (communications): 
     380      !    jj loop must be 1:jpj   if adv_x is called first 
     381      !                and 2:jpj-1 if adv_x is called second 
     382      jjmin = 2     - NINT(pcrh)   ! 1   or 2 
     383      jjmax = jpjm1 + NINT(pcrh)   ! jpj or jpj-1 
    276384      ! 
    277385      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    280388         ! 
    281389         ! Limitation of moments.                                            
    282          DO jj = 2, jpjm1 
     390         DO jj = jjmin, jjmax 
     391             
    283392            DO ji = 1, jpi 
     393 
     394               zpsm  = psm (ji,jj,jl) ! optimization 
     395               zps0  = ps0 (ji,jj,jl) 
     396               zpsx  = psx (ji,jj,jl) 
     397               zpsxx = psxx(ji,jj,jl) 
     398               zpsy  = psy (ji,jj,jl) 
     399               zpsyy = psyy(ji,jj,jl) 
     400               zpsxy = psxy(ji,jj,jl) 
     401 
    284402               !  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) ) 
     403               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     404               ! 
     405               zslpmax = MAX( 0._wp, zps0 ) 
    288406               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) )  ) 
     407               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     408               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
    292409               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    293410 
    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 
     411               zps0  = zslpmax   
     412               zpsx  = zs1new  * rswitch 
     413               zpsxx = zs2new  * rswitch 
     414               zpsy  = zpsy    * rswitch 
     415               zpsyy = zpsyy   * rswitch 
     416               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     417 
     418               !  Calculate fluxes and moments between boxes i<-->i+1               
     419               !                                !  Flux from i to i+1 WHEN u GT 0  
    306420               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) 
     421               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
    308422               zalfq        =  zalf * zalf 
    309423               zalf1        =  1.0 - zalf 
    310424               zalf1q       =  zalf1 * zalf1 
    311425               ! 
    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. 
     426               zfm (ji,jj)  =  zalf  *   zpsm  
     427               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     428               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     429               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     430               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     431               zfxy(ji,jj)  =  zalfq *   zpsxy 
     432               zfyy(ji,jj)  =  zalf  *   zpsyy 
     433 
     434               !                                !  Readjust moments remaining in the box. 
     435               zpsm  =  zpsm  - zfm(ji,jj) 
     436               zps0  =  zps0  - zf0(ji,jj) 
     437               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     438               zpsxx =  zalf1  * zalf1q * zpsxx 
     439               zpsy  =  zpsy  - zfy (ji,jj) 
     440               zpsyy =  zpsyy - zfyy(ji,jj) 
     441               zpsxy =  zalf1q * zpsxy 
     442               ! 
     443               psm (ji,jj,jl) = zpsm ! optimization 
     444               ps0 (ji,jj,jl) = zps0  
     445               psx (ji,jj,jl) = zpsx  
     446               psxx(ji,jj,jl) = zpsxx 
     447               psy (ji,jj,jl) = zpsy  
     448               psyy(ji,jj,jl) = zpsyy 
     449               psxy(ji,jj,jl) = zpsxy 
     450               ! 
     451            END DO 
     452 
    332453            DO ji = 1, fs_jpim1 
     454               !                                !  Flux from i+1 to i when u LT 0. 
    333455               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    334456               zalg  (ji,jj) = zalf 
    … …  
    348470               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    349471            END DO 
    350          END DO 
    351  
    352          DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.  
    353             DO ji = fs_2, fs_jpim1 
     472 
     473            DO ji = fs_2, fs_jpim1  
     474               ! 
     475               zpsm  = psm (ji,jj,jl) ! optimization 
     476               zps0  = ps0 (ji,jj,jl) 
     477               zpsx  = psx (ji,jj,jl) 
     478               zpsxx = psxx(ji,jj,jl) 
     479               zpsy  = psy (ji,jj,jl) 
     480               zpsyy = psyy(ji,jj,jl) 
     481               zpsxy = psxy(ji,jj,jl) 
     482               !                                !  Readjust moments remaining in the box. 
    354483               zbt  =       zbet(ji-1,jj) 
    355484               zbt1 = 1.0 - zbet(ji-1,jj) 
    356485               ! 
    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 
     486               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     487               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     488               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     489               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     490               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) ) 
     491               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     492               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     493 
     494               !   Put the temporary moments into appropriate neighboring boxes.     
     495               !                                !   Flux from i to i+1 IF u GT 0. 
     496               zbt   =       zbet(ji-1,jj) 
     497               zbt1  = 1.0 - zbet(ji-1,jj) 
     498               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     499               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     500               zalf1 = 1.0 - zalf 
     501               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     502               ! 
     503               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     504               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     505               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     506                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     507                  &            + zbt1 * zpsxx 
     508               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     509                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     510                  &            + zbt1 * zpsxy 
     511               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     512               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     513 
     514               !                                !  Flux from i+1 to i IF u LT 0. 
     515               zbt   =       zbet(ji,jj) 
     516               zbt1  = 1.0 - zbet(ji,jj) 
     517               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     518               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     519               zalf1 = 1.0 - zalf 
     520               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     521               ! 
     522               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     523               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     524               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     525                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     526                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     527               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     528                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     529               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     530               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     531               ! 
     532               psm (ji,jj,jl) = zpsm  ! optimization 
     533               ps0 (ji,jj,jl) = zps0  
     534               psx (ji,jj,jl) = zpsx  
     535               psxx(ji,jj,jl) = zpsxx 
     536               psy (ji,jj,jl) = zpsy  
     537               psyy(ji,jj,jl) = zpsyy 
     538               psxy(ji,jj,jl) = zpsxy 
     539            END DO 
     540             
    409541         END DO 
    410542 
    411543      END DO 
    412  
    413       !-- Lateral boundary conditions 
    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. ) 
    417544      ! 
    418545   END SUBROUTINE adv_x 
    … …  
    436563      !! 
    437564      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     565      INTEGER  ::   jimin, jimax                         ! dummy loop indices 
    438566      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    439567      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    440568      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     569      REAL(wp) ::   zpsm, zps0 
     570      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    441571      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    442572      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    443573      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    444574      !--------------------------------------------------------------------- 
     575      ! in order to avoid lbc_lnk (communications): 
     576      !    ji loop must be 1:jpi   if adv_y is called first 
     577      !                and 2:jpi-1 if adv_y is called second 
     578      jimin = 2     - NINT(pcrh)   ! 1   or 2 
     579      jimax = jpim1 + NINT(pcrh)   ! jpi or jpi-1 
    445580      ! 
    446581      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    450585         ! Limitation of moments. 
    451586         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) ) 
     587            DO ji = jimin, jimax 
     588               ! 
     589               zpsm  = psm (ji,jj,jl) ! optimization 
     590               zps0  = ps0 (ji,jj,jl) 
     591               zpsx  = psx (ji,jj,jl) 
     592               zpsxx = psxx(ji,jj,jl) 
     593               zpsy  = psy (ji,jj,jl) 
     594               zpsyy = psyy(ji,jj,jl) 
     595               zpsxy = psxy(ji,jj,jl) 
     596               ! 
     597               !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     598               zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     599               ! 
     600               zslpmax = MAX( 0._wp, zps0 ) 
    457601               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) )  ) 
     602               zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     603               zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    461604               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    462605               ! 
    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 
     606               zps0  = zslpmax   
     607               zpsx  = zpsx  * rswitch 
     608               zpsxx = zpsxx * rswitch 
     609               zpsy  = zs1new         * rswitch 
     610               zpsyy = zs2new         * rswitch 
     611               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
    471612  
    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 
     613               !  Calculate fluxes and moments between boxes j<-->j+1               
     614               !                                !  Flux from j to j+1 WHEN v GT 0    
    475615               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) 
     616               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    477617               zalfq        =  zalf * zalf 
    478618               zalf1        =  1.0 - zalf 
    479619               zalf1q       =  zalf1 * zalf1 
    480620               ! 
    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 
     621               zfm (ji,jj)  =  zalf  * zpsm 
     622               zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     623               zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     624               zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     625               zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     626               zfxy(ji,jj)  =  zalfq * zpsxy 
     627               zfxx(ji,jj)  =  zalf  * zpsxx 
     628               ! 
     629               !                                !  Readjust moments remaining in the box. 
     630               zpsm   =  zpsm  - zfm(ji,jj) 
     631               zps0   =  zps0  - zf0(ji,jj) 
     632               zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     633               zpsyy  =  zalf1 * zalf1q * zpsyy 
     634               zpsx   =  zpsx  - zfx(ji,jj) 
     635               zpsxx  =  zpsxx - zfxx(ji,jj) 
     636               zpsxy  =  zalf1q * zpsxy 
     637               ! 
     638               psm (ji,jj,jl) = zpsm ! optimization 
     639               ps0 (ji,jj,jl) = zps0  
     640               psx (ji,jj,jl) = zpsx  
     641               psxx(ji,jj,jl) = zpsxx 
     642               psy (ji,jj,jl) = zpsy  
     643               psyy(ji,jj,jl) = zpsyy 
     644               psxy(ji,jj,jl) = zpsxy 
     645            END DO 
     646         END DO 
     647         ! 
     648         DO jj = 1, jpjm1 
     649            DO ji = jimin, jimax 
     650               !                                !  Flux from j+1 to j when v LT 0. 
    502651               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    503652               zalg  (ji,jj) = zalf 
    … …  
    519668         END DO 
    520669 
    521          !  Readjust moments remaining in the box.  
    522670         DO jj = 2, jpjm1 
    523             DO ji = fs_2, fs_jpim1 
     671            DO ji = jimin, jimax 
     672               !                                !  Readjust moments remaining in the box. 
    524673               zbt  =         zbet(ji,jj-1) 
    525674               zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    526675               ! 
    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) ) 
     676               zpsm  = psm (ji,jj,jl) ! optimization 
     677               zps0  = ps0 (ji,jj,jl) 
     678               zpsx  = psx (ji,jj,jl) 
     679               zpsxx = psxx(ji,jj,jl) 
     680               zpsy  = psy (ji,jj,jl) 
     681               zpsyy = psyy(ji,jj,jl) 
     682               zpsxy = psxy(ji,jj,jl) 
     683               ! 
     684               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     685               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     686               zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     687               zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     688               zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     689               zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     690               zpsxy = zalg1q(ji,jj-1) * zpsxy 
     691 
     692               !   Put the temporary moments into appropriate neighboring boxes.     
     693               !                                !   Flux from j to j+1 IF v GT 0. 
     694               zbt   =       zbet(ji,jj-1) 
     695               zbt1  = 1.0 - zbet(ji,jj-1) 
     696               zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     697               zalf  = zbt * zfm(ji,jj-1) / zpsm  
     698               zalf1 = 1.0 - zalf 
     699               ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     700               ! 
     701               zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     702               zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     703                  &             + zbt1 * zpsy   
     704               zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     705                  &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     706                  &             + zbt1 * zpsyy 
     707               zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     708                  &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     709                  &             + zbt1 * zpsxy 
     710               zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     711               zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
     712 
     713               !                                !  Flux from j+1 to j IF v LT 0. 
     714               zbt   =       zbet(ji,jj) 
     715               zbt1  = 1.0 - zbet(ji,jj) 
     716               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     717               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     718               zalf1 = 1.0 - zalf 
     719               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     720               ! 
     721               zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     722               zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     723               zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     724                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     725                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     726               zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     727                  &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     728               zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     729               zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     730               ! 
     731               psm (ji,jj,jl) = zpsm ! optimization 
     732               ps0 (ji,jj,jl) = zps0  
     733               psx (ji,jj,jl) = zpsx  
     734               psxx(ji,jj,jl) = zpsxx 
     735               psy (ji,jj,jl) = zpsy  
     736               psyy(ji,jj,jl) = zpsyy 
     737               psxy(ji,jj,jl) = zpsxy 
    579738            END DO 
    580739         END DO 
    581740 
    582741      END DO 
    583  
    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. ) 
    588742      ! 
    589743   END SUBROUTINE adv_y 
     744 
     745 
     746   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     747      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
     748      !!------------------------------------------------------------------- 
     749      !!                  ***  ROUTINE Hbig  *** 
     750      !! 
     751      !! ** Purpose : Thickness correction in case advection scheme creates 
     752      !!              abnormally tick ice or snow 
     753      !! 
     754      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     755      !!                 (before advection) and reduce it by adapting ice concentration 
     756      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     757      !!                 (before advection) and reduce it by sending the excess in the ocean 
     758      !! 
     759      !! ** input   : Max thickness of the surrounding 9-points 
     760      !!------------------------------------------------------------------- 
     761      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     762      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     763      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     764      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     765      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
     766      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     767      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     768      ! 
     769      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     770      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
     771      !!------------------------------------------------------------------- 
     772      ! 
     773      z1_dt = 1._wp / pdt 
     774      ! 
     775      DO jl = 1, jpl 
     776         DO jj = 1, jpj 
     777            DO ji = 1, jpi 
     778               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     779                  ! 
     780                  !                               ! -- check h_ip -- ! 
     781                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     782                  IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     783                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     784                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     785                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     786                     ENDIF 
     787                  ENDIF 
     788                  ! 
     789                  !                               ! -- check h_i -- ! 
     790                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     791                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     792                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     793                     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) 
     794                  ENDIF 
     795                  ! 
     796                  !                               ! -- check h_s -- ! 
     797                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     798                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     799                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     800                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     801                     ! 
     802                     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 
     803                     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 
     804                     ! 
     805                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     806                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     807                  ENDIF            
     808                  !                   
     809                  !                               ! -- check s_i -- ! 
     810                  ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     811                  zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     812                  IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     813                     zfra = psi_max(ji,jj,jl) / zsi 
     814                     sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     815                     psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     816                  ENDIF 
     817                  ! 
     818               ENDIF 
     819            END DO 
     820         END DO 
     821      END DO  
     822      ! 
     823      !                                           ! -- check e_i/v_i -- ! 
     824      DO jl = 1, jpl 
     825         DO jk = 1, nlay_i 
     826            DO jj = 1, jpj 
     827               DO ji = 1, jpi 
     828                  IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     829                     ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     830                     zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     831                     IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     832                        zfra = pei_max(ji,jj,jk,jl) / zei 
     833                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     834                        pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     835                     ENDIF 
     836                  ENDIF 
     837               END DO 
     838            END DO 
     839         END DO 
     840      END DO 
     841      !                                           ! -- check e_s/v_s -- ! 
     842      DO jl = 1, jpl 
     843         DO jk = 1, nlay_s 
     844            DO jj = 1, jpj 
     845               DO ji = 1, jpi 
     846                  IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     847                     ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     848                     zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     849                     IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     850                        zfra = pes_max(ji,jj,jk,jl) / zes 
     851                        hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     852                        pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     853                     ENDIF 
     854                  ENDIF 
     855               END DO 
     856            END DO 
     857         END DO 
     858      END DO 
     859      ! 
     860   END SUBROUTINE Hbig 
    590861 
    591862 
    … …  
    659930         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   & 
    660931         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   & 
    661          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    662          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     932         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
     933         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     934         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   & 
    663935         ! 
    664936         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
    … …  
    7551027            END DO 
    7561028            ! 
    757             IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction 
    758                CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
    759                CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  ) 
    760                CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap ) 
    761                CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap ) 
    762                CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap ) 
    763                !                                                     ! melt pond volume 
    764                CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  ) 
    765                CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  ) 
    766                CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp ) 
    767                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
    768                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     1029            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     1030               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     1031                  CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
     1032                  CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap  ) 
     1033                  CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap ) 
     1034                  CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap ) 
     1035                  CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap ) 
     1036                  !                                                     ! melt pond volume 
     1037                  CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp  ) 
     1038                  CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp  ) 
     1039                  CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp ) 
     1040                  CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
     1041                  CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     1042               ELSE 
     1043                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction 
     1044                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume 
     1045               ENDIF 
     1046                  ! 
     1047               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
     1048                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     1049                     CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl  ) 
     1050                     CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl  ) 
     1051                     CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl ) 
     1052                     CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl ) 
     1053                     CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl ) 
     1054                  ELSE 
     1055                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume 
     1056                  ENDIF 
     1057               ENDIF 
    7691058            ENDIF 
    7701059            ! 
    … …  
    7801069            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    7811070            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    782             IF( ln_pnd_H12 ) THEN 
    783                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    784                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     1071            IF( ln_pnd_LEV ) THEN 
     1072               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
     1073               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
     1074               IF ( ln_pnd_lids ) THEN 
     1075                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume 
     1076               ENDIF 
    7851077            ENDIF 
    7861078         ENDIF 
    … …  
    8451137         END DO 
    8461138         ! 
    847          IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction 
     1139         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
    8481140            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    8491141            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 
    … …  
    8571149            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    8581150            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 
     1151            ! 
     1152            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume 
     1153               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  ) 
     1154               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  ) 
     1155               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl ) 
     1156               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl ) 
     1157               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl ) 
     1158            ENDIF 
    8591159         ENDIF 
    8601160         ! 
    … …  
    8631163   END SUBROUTINE adv_pra_rst 
    8641164 
     1165   SUBROUTINE icemax3D( pice , pmax ) 
     1166      !!--------------------------------------------------------------------- 
     1167      !!                   ***  ROUTINE icemax3D ***                      
     1168      !! ** Purpose :  compute the max of the 9 points around 
     1169      !!---------------------------------------------------------------------- 
     1170      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input 
     1171      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output 
     1172      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1173      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1174      !!---------------------------------------------------------------------- 
     1175      DO jl = 1, jpl 
     1176         DO jj = 1, jpj 
     1177            DO ji = 2, jpim1 
     1178               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 
     1179            END DO 
     1180         END DO 
     1181         DO jj = 2, jpjm1 
     1182            DO ji = 2, jpim1 
     1183               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1184            END DO 
     1185         END DO 
     1186      END DO 
     1187   END SUBROUTINE icemax3D 
     1188 
     1189   SUBROUTINE icemax4D( pice , pmax ) 
     1190      !!--------------------------------------------------------------------- 
     1191      !!                   ***  ROUTINE icemax4D ***                      
     1192      !! ** Purpose :  compute the max of the 9 points around 
     1193      !!---------------------------------------------------------------------- 
     1194      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input 
     1195      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output 
     1196      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array 
     1197      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices 
     1198      !!---------------------------------------------------------------------- 
     1199      jlay = SIZE( pice , 3 )   ! size of input arrays 
     1200      DO jl = 1, jpl 
     1201         DO jk = 1, jlay 
     1202            DO jj = 1, jpj 
     1203               DO ji = 2, jpim1 
     1204                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 
     1205               END DO 
     1206            END DO 
     1207            DO jj = 2, jpjm1 
     1208               DO ji = 2, jpim1 
     1209                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 
     1210               END DO 
     1211            END DO 
     1212         END DO 
     1213      END DO 
     1214   END SUBROUTINE icemax4D 
     1215    
    8651216#else 
    8661217   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.