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 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_pra.F90 – NEMO

Ignore:
Timestamp:
2020-12-03T12:20:38+01:00 (3 years ago)
Author:
ayoung
Message:

Updated to trunk at 14020. Sette tests passed with change of results for configurations with non-linear ssh. Ticket #2506.

Location:
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG

    • Property svn:externals

      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13292        sette 
         10^/utils/CI/sette_wave@13990         sette 

  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_adv_pra.F90

    r13295 r14037  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume 
    4647 
    4748   !! * Substitutions 
    … …  
    5556 
    5657   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
    … …  
    8182      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    8283      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness 
    8385      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8486      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8587      ! 
    86       INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     88      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8789      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    88       REAL(wp) ::   zdt                     !   -      - 
     90      REAL(wp) ::   zdt, z1_dt              !   -      - 
    8991      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
    9092      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    9193      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
    92       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     95      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    9397      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9498      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    95       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl 
    96100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    97101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     102      !! diagnostics 
     103      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat       
    98104      !!---------------------------------------------------------------------- 
    99105      ! 
    100106      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    101107      ! 
    102       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
    103       DO jl = 1, jpl 
    104          DO_2D( 0, 0, 0, 0 ) 
    105             zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    106                &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    107                &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    108                &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    109             zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    110                &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    111                &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    112                &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    113             zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    114                &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    115                &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    116                &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    117          END_2D 
     108      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     109      ! thickness and salinity 
     110      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     111      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     112      END WHERE 
     113      CALL icemax3D( ph_i , zhi_max ) 
     114      CALL icemax3D( ph_s , zhs_max ) 
     115      CALL icemax3D( ph_ip, zhip_max) 
     116      CALL icemax3D( zs_i , zsi_max ) 
     117      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
     118      ! 
     119      ! enthalpies 
     120      DO jk = 1, nlay_i 
     121         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     122         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     123         END WHERE 
    118124      END DO 
    119       CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 
     125      DO jk = 1, nlay_s 
     126         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     127         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     128         END WHERE 
     129      END DO    
     130      CALL icemax4D( ze_i , zei_max ) 
     131      CALL icemax4D( ze_s , zes_max ) 
     132      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 
     133      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 
     134      ! 
    120135      ! 
    121136      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
    … …  
    132147      ENDIF 
    133148      zdt = rDt_ice / REAL(icycle) 
     149      z1_dt = 1._wp / zdt 
    134150       
    135151      ! --- transport --- ! 
    … …  
    138154 
    139155      DO jt = 1, icycle 
     156 
     157         ! diagnostics 
     158         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     159            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
     160         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
     161         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     162            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 
    140163 
    141164         ! record at_i before advection (for open water) 
    … …  
    156179               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    157180            END DO 
    158             IF ( ln_pnd_H12 ) THEN 
    159                z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    160                z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     181            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
     182               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
     183               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
     184               IF ( ln_pnd_lids ) THEN 
     185                  z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume 
     186               ENDIF 
    161187            ENDIF 
    162188         END DO 
    … …  
    189215            END DO 
    190216            ! 
    191             IF ( ln_pnd_H12 ) THEN 
     217            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    192218               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    193219               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    194220               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    195221               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     222               IF ( ln_pnd_lids ) THEN 
     223                  CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     224                  CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     225               ENDIF 
    196226            ENDIF 
    197227            !                                                               !--------------------------------------------! 
    … …  
    220250                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    221251            END DO 
    222             IF ( ln_pnd_H12 ) THEN 
     252            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    223253               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    224254               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    225255               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    226256               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     257               IF ( ln_pnd_lids ) THEN 
     258                  CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     259                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     260               ENDIF 
    227261            ENDIF 
    228262            ! 
     263         ENDIF 
     264          
     265         ! --- Lateral boundary conditions --- ! 
     266         !     caution: for gradients (sx and sy) the sign changes 
     267         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume 
     268            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  & 
     269            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume 
     270            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  ) 
     271         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity 
     272            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  & 
     273            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration 
     274            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  ) 
     275         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age 
     276            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  ) 
     277         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy 
     278            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  )  
     279         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
     280            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
     281         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
     282            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
     283               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
     284               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume 
     285               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  )  
     286            IF ( ln_pnd_lids ) THEN 
     287               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume 
     288                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  )  
     289            ENDIF 
    229290         ENDIF 
    230291 
    … …  
    242303               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    243304            END DO 
    244             IF ( ln_pnd_H12 ) THEN 
     305            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    245306               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    246307               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     308               IF ( ln_pnd_lids ) THEN 
     309                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     310               ENDIF 
    247311            ENDIF 
    248312         END DO 
    … …  
    256320         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp ) 
    257321         ! 
     322         ! --- diagnostics --- ! 
     323         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     324            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
     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         ! 
    258332         ! --- Ensure non-negative fields --- ! 
    259333         !     Remove negative values (conservation is ensured) 
    260334         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    261          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     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 ) 
    262336         ! 
    263337         ! --- Make sure ice thickness is not too big --- ! 
    264338         !     (because ice thickness can be too large where ice concentration is very small) 
    265          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     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 ) 
    266341         ! 
    267342         ! --- Ensure snow load is not too big --- ! 
    … …  
    292367      !!  
    293368      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     369      INTEGER  ::   jj0                                  ! dummy loop indices 
    294370      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars 
    295371      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      - 
    296372      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      - 
     373      REAL(wp) ::   zpsm, zps0 
     374      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    297375      REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace 
    298376      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    299377      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    300378      !----------------------------------------------------------------------- 
     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      jj0 = NINT(pcrh) 
    301383      ! 
    302384      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    305387         ! 
    306388         ! Limitation of moments.                                            
    307          DO_2D( 0, 0, 1, 1 ) 
    308             !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    309             psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 
    310             ! 
    311             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
    312             zs1max  = 1.5 * zslpmax 
    313             zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 
    314             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    315                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  ) 
    316             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    317  
    318             ps0 (ji,jj,jl) = zslpmax   
    319             psx (ji,jj,jl) = zs1new         * rswitch 
    320             psxx(ji,jj,jl) = zs2new         * rswitch 
    321             psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 
    322             psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 
    323             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    324          END_2D 
    325  
    326          !  Calculate fluxes and moments between boxes i<-->i+1               
    327          DO_2D( 0, 0, 1, 1 ) 
    328             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    329             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    330             zalfq        =  zalf * zalf 
    331             zalf1        =  1.0 - zalf 
    332             zalf1q       =  zalf1 * zalf1 
    333             ! 
    334             zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl) 
    335             zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 
    336             zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 
    337             zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq 
    338             zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    339             zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl) 
    340             zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl) 
    341  
    342             !  Readjust moments remaining in the box. 
    343             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    344             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    345             psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 
    346             psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl) 
    347             psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj) 
    348             psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj) 
    349             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
    350          END_2D 
    351  
    352          DO_2D( 0, 0, 1, 0 ) 
    353             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    354             zalg  (ji,jj) = zalf 
    355             zalfq         = zalf * zalf 
    356             zalf1         = 1.0 - zalf 
    357             zalg1 (ji,jj) = zalf1 
    358             zalf1q        = zalf1 * zalf1 
    359             zalg1q(ji,jj) = zalf1q 
    360             ! 
    361             zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
    362             zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
    363                &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
    364             zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
    365             zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
    366             zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
    367             zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
    368             zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
    369          END_2D 
    370  
    371          DO_2D( 0, 0, 0, 0 ) 
    372             zbt  =       zbet(ji-1,jj) 
    373             zbt1 = 1.0 - zbet(ji-1,jj) 
    374             ! 
    375             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 
    376             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 
    377             psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 
    378             psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 
    379             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 
    380             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 
    381             psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 
    382          END_2D 
    383  
    384          !   Put the temporary moments into appropriate neighboring boxes.     
    385          DO_2D( 0, 0, 0, 0 ) 
    386             zbt  =       zbet(ji-1,jj) 
    387             zbt1 = 1.0 - zbet(ji-1,jj) 
    388             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 
    389             zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 
    390             zalf1         = 1.0 - zalf 
    391             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 
    392             ! 
    393             ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 
    394             psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 
    395             psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             & 
    396                &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) & 
    397                &            + zbt1 * psxx(ji,jj,jl) 
    398             psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             & 
    399                &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   & 
    400                &            + zbt1 * psxy(ji,jj,jl) 
    401             psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 
    402             psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 
    403          END_2D 
    404  
    405          DO_2D( 0, 0, 0, 0 ) 
    406             zbt  =       zbet(ji,jj) 
    407             zbt1 = 1.0 - zbet(ji,jj) 
    408             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    409             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    410             zalf1         = 1.0 - zalf 
    411             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    412             ! 
    413             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 
    414             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 
    415             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 
    416                &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    & 
    417                &                                           + ( zalf1 - zalf ) * ztemp ) ) 
    418             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    419                &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 
    420             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 
    421             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 
    422          END_2D 
    423  
     389         DO jj = Njs0 - jj0, Nje0 + jj0 
     390             
     391            DO ji = Nis0 - 1, Nie0 + 1 
     392 
     393               zpsm  = psm (ji,jj,jl) ! optimization 
     394               zps0  = ps0 (ji,jj,jl) 
     395               zpsx  = psx (ji,jj,jl) 
     396               zpsxx = psxx(ji,jj,jl) 
     397               zpsy  = psy (ji,jj,jl) 
     398               zpsyy = psyy(ji,jj,jl) 
     399               zpsxy = psxy(ji,jj,jl) 
     400 
     401               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
     402               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 
     403               ! 
     404               zslpmax = MAX( 0._wp, zps0 ) 
     405               zs1max  = 1.5 * zslpmax 
     406               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) ) 
     407               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 
     408               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
     409 
     410               zps0  = zslpmax   
     411               zpsx  = zs1new  * rswitch 
     412               zpsxx = zs2new  * rswitch 
     413               zpsy  = zpsy    * rswitch 
     414               zpsyy = zpsyy   * rswitch 
     415               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     416 
     417               !  Calculate fluxes and moments between boxes i<-->i+1               
     418               !                                !  Flux from i to i+1 WHEN u GT 0  
     419               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
     420               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 
     421               zalfq        =  zalf * zalf 
     422               zalf1        =  1.0 - zalf 
     423               zalf1q       =  zalf1 * zalf1 
     424               ! 
     425               zfm (ji,jj)  =  zalf  *   zpsm  
     426               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 
     427               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx ) 
     428               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq 
     429               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy ) 
     430               zfxy(ji,jj)  =  zalfq *   zpsxy 
     431               zfyy(ji,jj)  =  zalf  *   zpsyy 
     432 
     433               !                                !  Readjust moments remaining in the box. 
     434               zpsm  =  zpsm  - zfm(ji,jj) 
     435               zps0  =  zps0  - zf0(ji,jj) 
     436               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 
     437               zpsxx =  zalf1  * zalf1q * zpsxx 
     438               zpsy  =  zpsy  - zfy (ji,jj) 
     439               zpsyy =  zpsyy - zfyy(ji,jj) 
     440               zpsxy =  zalf1q * zpsxy 
     441               ! 
     442               psm (ji,jj,jl) = zpsm ! optimization 
     443               ps0 (ji,jj,jl) = zps0  
     444               psx (ji,jj,jl) = zpsx  
     445               psxx(ji,jj,jl) = zpsxx 
     446               psy (ji,jj,jl) = zpsy  
     447               psyy(ji,jj,jl) = zpsyy 
     448               psxy(ji,jj,jl) = zpsxy 
     449               ! 
     450            END DO 
     451             
     452            DO ji = Nis0 - 1, Nie0 
     453               !                                !  Flux from i+1 to i when u LT 0. 
     454               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
     455               zalg  (ji,jj) = zalf 
     456               zalfq         = zalf * zalf 
     457               zalf1         = 1.0 - zalf 
     458               zalg1 (ji,jj) = zalf1 
     459               zalf1q        = zalf1 * zalf1 
     460               zalg1q(ji,jj) = zalf1q 
     461               ! 
     462               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl) 
     463               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) & 
     464                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 
     465               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 
     466               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq 
     467               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 
     468               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl) 
     469               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl) 
     470            END DO 
     471 
     472            DO ji = Nis0, Nie0 
     473               ! 
     474               zpsm  = psm (ji,jj,jl) ! optimization 
     475               zps0  = ps0 (ji,jj,jl) 
     476               zpsx  = psx (ji,jj,jl) 
     477               zpsxx = psxx(ji,jj,jl) 
     478               zpsy  = psy (ji,jj,jl) 
     479               zpsyy = psyy(ji,jj,jl) 
     480               zpsxy = psxy(ji,jj,jl) 
     481               !                                !  Readjust moments remaining in the box. 
     482               zbt  =       zbet(ji-1,jj) 
     483               zbt1 = 1.0 - zbet(ji-1,jj) 
     484               ! 
     485               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 
     486               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 
     487               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 
     488               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 
     489               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) ) 
     490               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 
     491               zpsxy = zalg1q(ji-1,jj) * zpsxy 
     492 
     493               !   Put the temporary moments into appropriate neighboring boxes.     
     494               !                                !   Flux from i to i+1 IF u GT 0. 
     495               zbt   =       zbet(ji-1,jj) 
     496               zbt1  = 1.0 - zbet(ji-1,jj) 
     497               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 
     498               zalf  = zbt * zfm(ji-1,jj) / zpsm 
     499               zalf1 = 1.0 - zalf 
     500               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 
     501               ! 
     502               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 
     503               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 
     504               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            & 
     505                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 
     506                  &            + zbt1 * zpsxx 
     507               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            & 
     508                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  & 
     509                  &            + zbt1 * zpsxy 
     510               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy  
     511               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 
     512 
     513               !                                !  Flux from i+1 to i IF u LT 0. 
     514               zbt   =       zbet(ji,jj) 
     515               zbt1  = 1.0 - zbet(ji,jj) 
     516               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     517               zalf  = zbt1 * zfm(ji,jj) / zpsm 
     518               zalf1 = 1.0 - zalf 
     519               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     520               ! 
     521               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) ) 
     522               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 
     523               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 
     524                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    & 
     525                  &                         + ( zalf1 - zalf ) * ztemp ) ) 
     526               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     527                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 
     528               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) ) 
     529               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 
     530               ! 
     531               psm (ji,jj,jl) = zpsm  ! optimization 
     532               ps0 (ji,jj,jl) = zps0  
     533               psx (ji,jj,jl) = zpsx  
     534               psxx(ji,jj,jl) = zpsxx 
     535               psy (ji,jj,jl) = zpsy  
     536               psyy(ji,jj,jl) = zpsyy 
     537               psxy(ji,jj,jl) = zpsxy 
     538            END DO 
     539            ! 
     540         END DO 
     541         ! 
    424542      END DO 
    425  
    426       !-- Lateral boundary conditions 
    427       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    428          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    429          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    430       ! 
     543      !       
    431544   END SUBROUTINE adv_x 
    432545 
    … …  
    449562      !! 
    450563      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices 
     564      INTEGER  ::   ji0                                  ! dummy loop indices 
    451565      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars 
    452566      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         - 
    453567      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         - 
     568      REAL(wp) ::   zpsm, zps0 
     569      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy 
    454570      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace 
    455571      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      - 
    456572      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      - 
    457573      !--------------------------------------------------------------------- 
     574      ! in order to avoid lbc_lnk (communications): 
     575      !    ji loop must be 1:jpi   if adv_y is called first 
     576      !                and 2:jpi-1 if adv_y is called second 
     577      ji0 = NINT(pcrh) 
    458578      ! 
    459579      jcat = SIZE( ps0 , 3 )   ! size of input arrays 
    … …  
    462582         ! 
    463583         ! Limitation of moments. 
    464          DO_2D( 1, 1, 0, 0 ) 
    465             !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    466             psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  ) 
    467             ! 
    468             zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 
     584         DO_2D( 1, 1, ji0, ji0 ) 
     585            ! 
     586            zpsm  = psm (ji,jj,jl) ! optimization 
     587            zps0  = ps0 (ji,jj,jl) 
     588            zpsx  = psx (ji,jj,jl) 
     589            zpsxx = psxx(ji,jj,jl) 
     590            zpsy  = psy (ji,jj,jl) 
     591            zpsyy = psyy(ji,jj,jl) 
     592            zpsxy = psxy(ji,jj,jl) 
     593            ! 
     594            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 
     595            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  ) 
     596            ! 
     597            zslpmax = MAX( 0._wp, zps0 ) 
    469598            zs1max  = 1.5 * zslpmax 
    470             zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 
    471             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    472                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  ) 
     599            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) ) 
     600            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 
    473601            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask 
    474602            ! 
    475             ps0 (ji,jj,jl) = zslpmax   
    476             psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 
    477             psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 
    478             psy (ji,jj,jl) = zs1new         * rswitch 
    479             psyy(ji,jj,jl) = zs2new         * rswitch 
    480             psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 
    481          END_2D 
    482   
    483          !  Calculate fluxes and moments between boxes j<-->j+1               
    484          DO_2D( 1, 1, 0, 0 ) 
     603            zps0  = zslpmax   
     604            zpsx  = zpsx  * rswitch 
     605            zpsxx = zpsxx * rswitch 
     606            zpsy  = zs1new         * rswitch 
     607            zpsyy = zs2new         * rswitch 
     608            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 
     609 
     610            !  Calculate fluxes and moments between boxes j<-->j+1               
     611            !                                !  Flux from j to j+1 WHEN v GT 0    
    485612            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    486             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     613            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 
    487614            zalfq        =  zalf * zalf 
    488615            zalf1        =  1.0 - zalf 
    489616            zalf1q       =  zalf1 * zalf1 
    490617            ! 
    491             zfm (ji,jj)  =  zalf  * psm(ji,jj,jl) 
    492             zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) )  
    493             zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 
    494             zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl) 
    495             zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 
    496             zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl) 
    497             zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl) 
    498             ! 
    499             !  Readjust moments remaining in the box. 
    500             psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj) 
    501             ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj) 
    502             psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 
    503             psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl) 
    504             psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj) 
    505             psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj) 
    506             psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl) 
     618            zfm (ji,jj)  =  zalf  * zpsm 
     619            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) )  
     620            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy ) 
     621            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy 
     622            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy ) 
     623            zfxy(ji,jj)  =  zalfq * zpsxy 
     624            zfxx(ji,jj)  =  zalf  * zpsxx 
     625            ! 
     626            !                                !  Readjust moments remaining in the box. 
     627            zpsm   =  zpsm  - zfm(ji,jj) 
     628            zps0   =  zps0  - zf0(ji,jj) 
     629            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 
     630            zpsyy  =  zalf1 * zalf1q * zpsyy 
     631            zpsx   =  zpsx  - zfx(ji,jj) 
     632            zpsxx  =  zpsxx - zfxx(ji,jj) 
     633            zpsxy  =  zalf1q * zpsxy 
     634            ! 
     635            psm (ji,jj,jl) = zpsm ! optimization 
     636            ps0 (ji,jj,jl) = zps0  
     637            psx (ji,jj,jl) = zpsx  
     638            psxx(ji,jj,jl) = zpsxx 
     639            psy (ji,jj,jl) = zpsy  
     640            psyy(ji,jj,jl) = zpsyy 
     641            psxy(ji,jj,jl) = zpsxy 
    507642         END_2D 
    508643         ! 
    509          DO_2D( 1, 0, 0, 0 ) 
     644         DO_2D( 1, 0, ji0, ji0 ) 
     645            !                                !  Flux from j+1 to j when v LT 0. 
    510646            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    511647            zalg  (ji,jj) = zalf 
    … …  
    526662         END_2D 
    527663 
    528          !  Readjust moments remaining in the box.  
    529          DO_2D( 0, 0, 0, 0 ) 
     664         DO_2D( 0, 0, ji0, ji0 ) 
     665            !                                !  Readjust moments remaining in the box. 
    530666            zbt  =         zbet(ji,jj-1) 
    531667            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    532668            ! 
    533             psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 
    534             ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 
    535             psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 
    536             psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 
    537             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 
    538             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 
    539             psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 
     669            zpsm  = psm (ji,jj,jl) ! optimization 
     670            zps0  = ps0 (ji,jj,jl) 
     671            zpsx  = psx (ji,jj,jl) 
     672            zpsxx = psxx(ji,jj,jl) 
     673            zpsy  = psy (ji,jj,jl) 
     674            zpsyy = psyy(ji,jj,jl) 
     675            zpsxy = psxy(ji,jj,jl) 
     676            ! 
     677            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 
     678            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 
     679            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 
     680            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 
     681            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) ) 
     682            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 
     683            zpsxy = zalg1q(ji,jj-1) * zpsxy 
     684 
     685            !   Put the temporary moments into appropriate neighboring boxes.     
     686            !                                !   Flux from j to j+1 IF v GT 0. 
     687            zbt   =       zbet(ji,jj-1) 
     688            zbt1  = 1.0 - zbet(ji,jj-1) 
     689            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm  
     690            zalf  = zbt * zfm(ji,jj-1) / zpsm  
     691            zalf1 = 1.0 - zalf 
     692            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 
     693            ! 
     694            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 
     695            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  & 
     696               &             + zbt1 * zpsy   
     697            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           & 
     698               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
     699               &             + zbt1 * zpsyy 
     700            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             & 
     701               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  & 
     702               &             + zbt1 * zpsxy 
     703            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx  
     704            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 
     705 
     706            !                                !  Flux from j+1 to j IF v LT 0. 
     707            zbt   =       zbet(ji,jj) 
     708            zbt1  = 1.0 - zbet(ji,jj) 
     709            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 
     710            zalf  = zbt1 * zfm(ji,jj) / zpsm 
     711            zalf1 = 1.0 - zalf 
     712            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 
     713            ! 
     714            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) ) 
     715            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 
     716            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 
     717               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     & 
     718               &                         + ( zalf1 - zalf ) * ztemp ) ) 
     719            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  & 
     720               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 
     721            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) ) 
     722            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 
     723            ! 
     724            psm (ji,jj,jl) = zpsm ! optimization 
     725            ps0 (ji,jj,jl) = zps0  
     726            psx (ji,jj,jl) = zpsx  
     727            psxx(ji,jj,jl) = zpsxx 
     728            psy (ji,jj,jl) = zpsy  
     729            psyy(ji,jj,jl) = zpsyy 
     730            psxy(ji,jj,jl) = zpsxy 
    540731         END_2D 
    541  
    542          !   Put the temporary moments into appropriate neighboring boxes.     
    543          DO_2D( 0, 0, 0, 0 ) 
    544             zbt  =       zbet(ji,jj-1) 
    545             zbt1 = 1.0 - zbet(ji,jj-1) 
    546             psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)  
    547             zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)  
    548             zalf1         = 1.0 - zalf 
    549             ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 
    550             ! 
    551             ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 
    552             psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  & 
    553                &             + zbt1 * psy(ji,jj,jl)   
    554             psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           & 
    555                &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &  
    556                &             + zbt1 * psyy(ji,jj,jl) 
    557             psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            & 
    558                &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  & 
    559                &             + zbt1 * psxy(ji,jj,jl) 
    560             psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 
    561             psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 
    562          END_2D 
    563  
    564          DO_2D( 0, 0, 0, 0 ) 
    565             zbt  =       zbet(ji,jj) 
    566             zbt1 = 1.0 - zbet(ji,jj) 
    567             psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 
    568             zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 
    569             zalf1         = 1.0 - zalf 
    570             ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 
    571             ! 
    572             ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) ) 
    573             psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 
    574             psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 
    575                &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    & 
    576                &                                            + ( zalf1 - zalf ) * ztemp ) ) 
    577             psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  & 
    578                &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 
    579             psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 
    580             psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 
    581          END_2D 
    582  
     732         ! 
    583733      END DO 
    584  
    585       !-- Lateral boundary conditions 
    586       CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   & 
    587          &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes 
    588          &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp ) 
    589734      ! 
    590735   END SUBROUTINE adv_y 
    591736 
    592737 
    593    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     738   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     739      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    594740      !!------------------------------------------------------------------- 
    595741      !!                  ***  ROUTINE Hbig  *** 
    … …  
    605751      !! ** input   : Max thickness of the surrounding 9-points 
    606752      !!------------------------------------------------------------------- 
    607       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    608       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    609       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     753      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     754      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     755      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     756      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     757      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    610758      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    611       ! 
    612       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    613       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     759      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     760      ! 
     761      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     762      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    614763      !!------------------------------------------------------------------- 
    615764      ! 
    … …  
    617766      ! 
    618767      DO jl = 1, jpl 
    619  
    620768         DO_2D( 1, 1, 1, 1 ) 
    621769            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    … …  
    623771               !                               ! -- check h_ip -- ! 
    624772               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    625                IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     773               IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    626774                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    627775                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    … …  
    650798               ENDIF            
    651799               !                   
     800               !                               ! -- check s_i -- ! 
     801               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     802               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     803               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     804                  zfra = psi_max(ji,jj,jl) / zsi 
     805                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     806                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     807               ENDIF 
     808               ! 
    652809            ENDIF 
    653810         END_2D 
    654811      END DO  
     812      ! 
     813      !                                           ! -- check e_i/v_i -- ! 
     814      DO jl = 1, jpl 
     815         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     816            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     817               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     818               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     819               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     820                  zfra = pei_max(ji,jj,jk,jl) / zei 
     821                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     822                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     823               ENDIF 
     824            ENDIF 
     825         END_3D 
     826      END DO 
     827      !                                           ! -- check e_s/v_s -- ! 
     828      DO jl = 1, jpl 
     829         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     830            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     831               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     832               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     833               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     834                  zfra = pes_max(ji,jj,jk,jl) / zes 
     835                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     836                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     837               ENDIF 
     838            ENDIF 
     839         END_3D 
     840      END DO 
    655841      ! 
    656842   END SUBROUTINE Hbig 
    … …  
    724910         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   & 
    725911         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   & 
    726          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    727          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     912         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
     913         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     914         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   & 
    728915         ! 
    729916         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
    … …  
    772959            ! 
    773960            !                                                        ! ice thickness 
    774             CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice ) 
    775             CALL iom_get( numrir, jpdom_auto, 'syice' , syice ) 
     961            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 
     962            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 
    776963            CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 
    777964            CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 
    778965            CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 
    779966            !                                                        ! snow thickness 
    780             CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  ) 
    781             CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  ) 
     967            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp ) 
     968            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp ) 
    782969            CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  ) 
    783970            CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  ) 
    784971            CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  ) 
    785972            !                                                        ! ice concentration 
    786             CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    ) 
    787             CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    ) 
     973            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp ) 
     974            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp ) 
    788975            CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   ) 
    789976            CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   ) 
    790977            CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   ) 
    791978            !                                                        ! ice salinity 
    792             CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal ) 
    793             CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal ) 
     979            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 
     980            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 
    794981            CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 
    795982            CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 
    796983            CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 
    797984            !                                                        ! ice age 
    798             CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage ) 
    799             CALL iom_get( numrir, jpdom_auto, 'syage' , syage ) 
     985            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 
     986            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 
    800987            CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 
    801988            CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 
    … …  
    804991            DO jk = 1, nlay_s 
    805992               WRITE(zchar1,'(I2.2)') jk 
    806                znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
    807                znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
    808                znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
    809                znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
    810                znam = 'sxyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
     993               znam = 'sxc0'//'_l'//zchar1   
     994               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:) 
     995               znam = 'syc0'//'_l'//zchar1   
     996               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:) 
     997               znam = 'sxxc0'//'_l'//zchar1  
     998               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:) 
     999               znam = 'syyc0'//'_l'//zchar1  
     1000               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:) 
     1001               znam = 'sxyc0'//'_l'//zchar1  
     1002               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:) 
    8111003            END DO 
    8121004            !                                                        ! ice layers heat content 
    8131005            DO jk = 1, nlay_i 
    8141006               WRITE(zchar1,'(I2.2)') jk 
    815                znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
    816                znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
    817                znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
    818                znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
    819                znam = 'sxye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
    820             END DO 
    821             ! 
    822             IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction 
    823                CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  ) 
    824                CALL iom_get( numrir, jpdom_auto, 'syap' , syap  ) 
    825                CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
    826                CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
    827                CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
    828                !                                                     ! melt pond volume 
    829                CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  ) 
    830                CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  ) 
    831                CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
    832                CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
    833                CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
     1007               znam = 'sxe'//'_l'//zchar1    
     1008               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:) 
     1009               znam = 'sye'//'_l'//zchar1    
     1010               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:) 
     1011               znam = 'sxxe'//'_l'//zchar1   
     1012               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:) 
     1013               znam = 'syye'//'_l'//zchar1   
     1014               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:) 
     1015               znam = 'sxye'//'_l'//zchar1   
     1016               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:) 
     1017            END DO 
     1018            ! 
     1019            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction 
     1020               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     1021                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 
     1022                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 
     1023                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
     1024                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
     1025                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
     1026                  !                                                     ! melt pond volume 
     1027                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 
     1028                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 
     1029                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
     1030                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
     1031                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
     1032               ELSE 
     1033                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction 
     1034                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume 
     1035               ENDIF 
     1036                  ! 
     1037               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
     1038                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     1039                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 
     1040                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 
     1041                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 
     1042                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 
     1043                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl ) 
     1044                  ELSE 
     1045                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume 
     1046                  ENDIF 
     1047               ENDIF 
    8341048            ENDIF 
    8351049            ! 
    … …  
    8451059            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    8461060            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    847             IF( ln_pnd_H12 ) THEN 
    848                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    849                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     1061            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
     1062               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
     1063               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
     1064               IF ( ln_pnd_lids ) THEN 
     1065                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume 
     1066               ENDIF 
    8501067            ENDIF 
    8511068         ENDIF 
    … …  
    8621079         ! 
    8631080         !                                                           ! ice thickness 
    864          CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  ) 
    865          CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  ) 
    866          CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice ) 
    867          CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice ) 
    868          CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice ) 
     1081         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice) 
     1082         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice) 
     1083         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice) 
     1084         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice) 
     1085         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice) 
    8691086         !                                                           ! snow thickness 
    870          CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   ) 
    871          CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   ) 
    872          CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  ) 
    873          CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  ) 
    874          CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  ) 
     1087         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn ) 
     1088         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn ) 
     1089         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn ) 
     1090         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn ) 
     1091         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn ) 
    8751092         !                                                           ! ice concentration 
    876          CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    ) 
    877          CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    ) 
    878          CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   ) 
    879          CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   ) 
    880          CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   ) 
     1093         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa  ) 
     1094         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya  ) 
     1095         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa  ) 
     1096         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya  ) 
     1097         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya  ) 
    8811098         !                                                           ! ice salinity 
    882          CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  ) 
    883          CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  ) 
    884          CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal ) 
    885          CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal ) 
    886          CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal ) 
     1099         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal) 
     1100         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal) 
     1101         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal) 
     1102         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal) 
     1103         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal) 
    8871104         !                                                           ! ice age 
    888          CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  ) 
    889          CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  ) 
    890          CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage ) 
    891          CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 
    892          CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 
     1105         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage) 
     1106         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage) 
     1107         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage) 
     1108         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage) 
     1109         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage) 
    8931110         !                                                           ! snow layers heat content 
    8941111         DO jk = 1, nlay_s 
    8951112            WRITE(zchar1,'(I2.2)') jk 
    896             znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    897             znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    898             znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    899             znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    900             znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)  ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
     1113            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:) 
     1114            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1115            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:) 
     1116            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1117            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:) 
     1118            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1119            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:) 
     1120            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1121            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:) 
     1122            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
    9011123         END DO 
    9021124         !                                                           ! ice layers heat content 
    9031125         DO jk = 1, nlay_i 
    9041126            WRITE(zchar1,'(I2.2)') jk 
    905             znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    906             znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    907             znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    908             znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
    909             znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)   ;   CALL iom_rstput( iter, nitrst, numriw, znam , z3d ) 
     1127            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:) 
     1128            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1129            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:) 
     1130            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1131            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:) 
     1132            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1133            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:) 
     1134            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
     1135            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:) 
     1136            CALL iom_rstput( iter, nitrst, numriw, znam , z3d) 
    9101137         END DO 
    9111138         ! 
    912          IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction 
     1139         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction 
    9131140            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    9141141            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 
    … …  
    9221149            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    9231150            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 
    9241159         ENDIF 
    9251160         ! 
    … …  
    9271162      ! 
    9281163   END SUBROUTINE adv_pra_rst 
     1164 
     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 = Njs0-1, Nje0+1     
     1177            DO ji = Nis0, Nie0 
     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 = Njs0, Nje0     
     1182            DO ji = Nis0, Nie0 
     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 = Njs0-1, Nje0+1     
     1203               DO ji = Nis0, Nie0 
     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 = Njs0, Nje0     
     1208               DO ji = Nis0, Nie0 
     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 
    9291215 
    9301216#else 
Note: See TracChangeset for help on using the changeset viewer.