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 5709 for branches/UKMO/DEV_r5107_dynvor_updates/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 – NEMO

Ignore:
Timestamp:
2015-08-26T16:14:53+02:00 (9 years ago)
Author:
davestorkey
Message:

Updating UKMO/DEV_r5107_dynvor_updates branch to be relative to revision 5518 of the trunk
(= branching point for NEMO 3.6_STABLE).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/DEV_r5107_dynvor_updates/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5256 r5709  
    1717   USE dom_oce        ! ocean domain 
    1818   USE sbc_oce        ! ocean surface boundary condition 
    19    USE par_ice        ! ice parameter 
    2019   USE dom_ice        ! ice domain 
    2120   USE ice            ! ice variables 
    2221   USE limadv         ! ice advection 
    2322   USE limhdf         ! ice horizontal diffusion 
     23   USE limvar         !  
     24   ! 
    2425   USE in_out_manager ! I/O manager 
    2526   USE lbclnk         ! lateral boundary conditions -- MPP exchanges 
     
    2829   USE prtctl         ! Print control 
    2930   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    30    USE limvar          ! clem for ice thickness correction 
    31    USE timing          ! Timing 
     31   USE timing         ! Timing 
    3232   USE limcons        ! conservation tests 
     33   USE limctl         ! control prints 
    3334 
    3435   IMPLICIT NONE 
    3536   PRIVATE 
    3637 
    37    PUBLIC   lim_trp    ! called by ice_step 
     38   PUBLIC   lim_trp    ! called by sbcice_lim 
     39 
     40   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2   
    3841 
    3942   !! * Substitution 
     
    5861      !! ** action : 
    5962      !!--------------------------------------------------------------------- 
    60       INTEGER, INTENT(in) ::   kt   ! number of iteration 
     63      INTEGER, INTENT(in) ::   kt           ! number of iteration 
    6164      ! 
    62       INTEGER  ::   ji, jj, jk, jl, jn      ! dummy loop indices 
     65      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    6366      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    6467      REAL(wp) ::   zcfl , zusnit           !   -      - 
     68      CHARACTER(len=80) ::   cltmp 
    6569      ! 
    66       REAL(wp), POINTER, DIMENSION(:,:)      ::   zui_u, zvi_v, zsm, zs0at, zs0ow 
    67       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 
    68       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
    69       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold   ! old ice volume... 
    70       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zaiold, zhimax   ! old ice concentration and thickness 
    71       REAL(wp), POINTER, DIMENSION(:,:)      ::   zeiold, zesold   ! old enthalpies 
    72       REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 
    73       ! 
    74       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     70      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm 
     71      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
     72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw 
     73      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei 
     74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume... 
     75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness 
     76      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies 
     77      REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei 
     78      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
    7579      !!--------------------------------------------------------------------- 
    7680      IF( nn_timing == 1 )  CALL timing_start('limtrp') 
    7781 
    78       CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 
    79       CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 
    80       CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 
    81  
    82       CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem 
     82      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold ) 
     83      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     84      CALL wrk_alloc( jpi,jpj,1,          z0opw ) 
     85      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     86      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold ) 
    8387 
    8488      IF( numit == nstart .AND. lwp ) THEN 
     
    8892         ENDIF 
    8993         WRITE(numout,*) '~~~~~~~~~~~~' 
     94         ncfl = 0                ! nb of time step with CFL > 1/2 
    9095      ENDIF 
     96 
     97      zsm(:,:) = e12t(:,:) 
    9198       
    92       zsm(:,:) = area(:,:) 
    93  
    9499      !                             !-------------------------------------! 
    95100      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   ! 
     
    97102 
    98103         ! conservation test 
    99          IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    100  
    101          ! mass and salt flux init (clem) 
     104         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     105 
     106         ! mass and salt flux init 
    102107         zviold(:,:,:)  = v_i(:,:,:) 
    103          zeiold(:,:)  = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
    104          zesold(:,:)  = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
    105  
    106          !--- Thickness correction init. (clem) ------------------------------- 
    107          CALL lim_var_glo2eqv 
    108          zaiold(:,:,:) = a_i(:,:,:) 
     108         zvsold(:,:,:)  = v_s(:,:,:) 
     109         zsmvold(:,:,:) = smv_i(:,:,:) 
     110         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
     111         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
     112 
     113         !--- Thickness correction init. ------------------------------- 
     114         zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     115         DO jl = 1, jpl 
     116            DO jj = 1, jpj 
     117               DO ji = 1, jpi 
     118                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     119                  ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     120                  ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     121               END DO 
     122            END DO 
     123         END DO 
    109124         !--------------------------------------------------------------------- 
    110          ! Record max of the surrounding ice thicknesses for correction in limupdate 
     125         ! Record max of the surrounding ice thicknesses for correction 
    111126         ! in case advection creates ice too thick. 
    112127         !--------------------------------------------------------------------- 
    113          zhimax(:,:,:) = ht_i(:,:,:) 
     128         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 
    114129         DO jl = 1, jpl 
    115130            DO jj = 2, jpjm1 
    116131               DO ji = 2, jpim1 
    117                   zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 
    118                   !zhimax(ji,jj,jl) = ( ht_i(ji  ,jj  ,jl) * tmask(ji,  jj  ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) & 
    119                   !     &             + ht_i(ji-1,jj  ,jl) * tmask(ji-1,jj  ,1) + ht_i(ji  ,jj-1,jl) * tmask(ji  ,jj-1,1) & 
    120                   !     &             + ht_i(ji+1,jj  ,jl) * tmask(ji+1,jj  ,1) + ht_i(ji  ,jj+1,jl) * tmask(ji  ,jj+1,1) & 
    121                   !     &             + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) ) 
     132                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 
    122133               END DO 
    123134            END DO 
     
    125136         END DO 
    126137          
     138         !=============================! 
     139         !==      Prather scheme     ==! 
     140         !=============================! 
     141 
     142         ! If ice drift field is too fast, use an appropriate time step for advection.          
     143         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability 
     144         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
     145         IF(lk_mpp )   CALL mpp_max( zcfl ) 
     146 
     147         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
     148         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     149         ENDIF 
     150 
     151         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1 
     152!!         IF( lwp ) THEN 
     153!!            IF( ncfl > 0 ) THEN    
     154!!               WRITE(cltmp,'(i6.1)') ncfl 
     155!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
     156!!            ELSE 
     157!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 
     158!!            ENDIF 
     159!!         ENDIF 
     160 
    127161         !------------------------- 
    128162         ! transported fields                                         
    129163         !------------------------- 
    130          ! Snow vol, ice vol, salt and age contents, area 
    131          zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area  
    132          DO jl = 1, jpl 
    133             zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume 
    134             zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume 
    135             zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area 
    136             zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content 
    137             zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content 
    138             zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content 
    139             zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content 
    140          END DO 
    141  
    142          !-------------------------- 
    143          ! Advection of Ice fields  (Prather scheme)                                             
    144          !-------------------------- 
    145          ! If ice drift field is too fast, use an appropriate time step for advection.          
    146          ! CFL test for stability 
    147          zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 
    148          zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 
    149          IF(lk_mpp )   CALL mpp_max( zcfl ) 
    150 !!gm more readability: 
    151 !         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
    152 !         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
    153 !         ENDIF 
    154 !!gm end 
    155          initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 
    156          zusnit = 1.0 / REAL( initad )  
    157          IF( zcfl > 0.5 .AND. lwp )   & 
    158             WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
    159                &                        ': the ice time stepping is split in two' 
     164         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area  
     165         DO jl = 1, jpl 
     166            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
     167            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
     168            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
     169            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
     170            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
     171            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
     172            DO jk = 1, nlay_i 
     173               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
     174            END DO 
     175         END DO 
     176 
    160177 
    161178         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    162             DO jn = 1,initad 
    163                CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
    164                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    165                CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   & 
    166                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     179            DO jt = 1, initad 
     180               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     181                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     182               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   & 
     183                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    167184               DO jl = 1, jpl 
    168                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     185                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    169186                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    170                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     187                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   & 
    171188                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    172                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     189                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    173190                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    174                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     191                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    175192                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    176                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     193                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    177194                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    178                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     195                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    179196                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    180                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---      
     197                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    181198                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    182                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     199                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   & 
    183200                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    184                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     201                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    185202                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    186                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &  
     203                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    187204                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    188                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     205                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
    189206                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    190                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     207                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   & 
    191208                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    192                   DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    193                      CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     209                  DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
     210                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    194211                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    195212                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    196                      CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     213                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    197214                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    198215                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     
    201218            END DO 
    202219         ELSE 
    203             DO jn = 1, initad 
    204                CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
    205                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    206                CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   & 
    207                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     220            DO jt = 1, initad 
     221               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     222                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     223               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   & 
     224                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    208225               DO jl = 1, jpl 
    209                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     226                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    210227                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    211                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     228                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   & 
    212229                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    213                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     230                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    214231                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    215                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     232                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    216233                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    217                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     234                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    218235                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    219                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     236                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    220237                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    221  
    222                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     238                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    223239                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    224                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     240                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   & 
    225241                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    226                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     242                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    227243                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    228                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
     244                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    229245                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    230                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     246                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    231247                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    232                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     248                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   & 
    233249                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    234250                  DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    235                      CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     251                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    236252                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    237253                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    238                      CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     254                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    239255                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    240256                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
     
    247263         ! Recover the properties from their contents 
    248264         !------------------------------------------- 
    249          zs0ow(:,:) = zs0ow(:,:) / area(:,:) 
    250          DO jl = 1, jpl 
    251             zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) 
    252             zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:) 
    253             zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:) 
    254             zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 
    255             zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:) 
    256             ! 
     265         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 
     266         DO jl = 1, jpl 
     267            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:) 
     268            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:) 
     269            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:) 
     270            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:) 
     271            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:) 
     272            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 
     273            DO jk = 1, nlay_i 
     274               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 
     275            END DO 
     276         END DO 
     277 
     278         at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
     279         DO jl = 2, jpl 
     280            at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    257281         END DO 
    258282 
    259283         !------------------------------------------------------------------------------! 
    260          ! 4) Diffusion of Ice fields                   
     284         ! Diffusion of Ice fields                   
    261285         !------------------------------------------------------------------------------! 
    262286 
     287         ! 
    263288         !-------------------------------- 
    264289         !  diffusion of open water area 
    265290         !-------------------------------- 
    266          zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction 
    267          DO jl = 2, jpl 
    268             zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl) 
    269          END DO 
    270          ! 
    271291         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    272292         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
    273293            DO ji = 1 , fs_jpim1   ! vector opt. 
    274                pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji  ,jj) ) ) )   & 
    275                   &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 
    276                pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj  ) ) ) )   & 
    277                   &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 
     294               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   & 
     295                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 
     296               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   & 
     297                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 
    278298            END DO 
    279299         END DO 
    280300         ! 
    281          CALL lim_hdf( zs0ow (:,:) )   ! Diffusion 
     301         CALL lim_hdf( ato_i (:,:) ) 
    282302 
    283303         !------------------------------------ 
     
    288308            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
    289309               DO ji = 1 , fs_jpim1   ! vector opt. 
    290                   pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji  ,jj,jl) ) ) )   & 
    291                      &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
    292                   pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj  ,jl) ) ) )   & 
    293                      &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
    294                END DO 
    295             END DO 
    296  
    297             CALL lim_hdf( zs0ice (:,:,jl) ) 
    298             CALL lim_hdf( zs0sn  (:,:,jl) ) 
    299             CALL lim_hdf( zs0sm  (:,:,jl) ) 
    300             CALL lim_hdf( zs0oi  (:,:,jl) ) 
    301             CALL lim_hdf( zs0a   (:,:,jl) ) 
    302             CALL lim_hdf( zs0c0  (:,:,jl) ) 
     310                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   & 
     311                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
     312                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   & 
     313                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
     314               END DO 
     315            END DO 
     316 
     317            CALL lim_hdf( v_i  (:,:,  jl) ) 
     318            CALL lim_hdf( v_s  (:,:,  jl) ) 
     319            CALL lim_hdf( smv_i(:,:,  jl) ) 
     320            CALL lim_hdf( oa_i (:,:,  jl) ) 
     321            CALL lim_hdf( a_i  (:,:,  jl) ) 
     322            CALL lim_hdf( e_s  (:,:,1,jl) ) 
    303323            DO jk = 1, nlay_i 
    304                CALL lim_hdf( zs0e (:,:,jk,jl) ) 
     324               CALL lim_hdf( e_i(:,:,jk,jl) ) 
    305325            END DO 
    306326         END DO 
    307327 
    308328         !------------------------------------------------------------------------------! 
    309          ! 5) Update and limit ice properties after transport                            
     329         ! limit ice properties after transport                            
    310330         !------------------------------------------------------------------------------! 
    311  
    312          !-------------------------------------------------- 
    313          ! 5.1) Recover mean values over the grid squares. 
    314          !-------------------------------------------------- 
    315          zs0at(:,:) = 0._wp 
     331!!gm & cr   :  MAX should not be active if adv scheme is positive ! 
    316332         DO jl = 1, jpl 
    317333            DO jj = 1, jpj 
    318334               DO ji = 1, jpi 
    319                   zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 
    320                   zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 
    321                   zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 
    322                   zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 
    323                   zs0a  (ji,jj,jl) = MAX( 0._wp, zs0a  (ji,jj,jl) ) 
    324                   zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 
    325                   zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl) 
    326                END DO 
    327             END DO 
    328          END DO 
    329  
    330          !--------------------------------------------------------- 
    331          ! 5.2) Update and mask variables 
    332          !--------------------------------------------------------- 
    333          DO jl = 1, jpl           
    334             DO jj = 1, jpj 
    335                DO ji = 1, jpi 
    336                   rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 
    337  
    338                   zvi  = zs0ice(ji,jj,jl) 
    339                   zvs  = zs0sn (ji,jj,jl) 
    340                   zes  = zs0c0 (ji,jj,jl)       
    341                   zsmv = zs0sm (ji,jj,jl) 
    342                   ! 
    343                   ! Remove very small areas 
    344                   v_s(ji,jj,jl)   = rswitch * zs0sn (ji,jj,jl)  
    345                   v_i(ji,jj,jl)   = rswitch * zs0ice(ji,jj,jl) 
    346                   a_i(ji,jj,jl)   = rswitch * zs0a  (ji,jj,jl) 
    347                   e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl)       
    348                   ! Ice salinity and age 
    349                   IF(  num_sal == 2  ) THEN 
    350                      smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 
    351                   ENDIF 
    352                   oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 
    353  
    354                  ! Update fluxes 
    355                   wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice  
    356                   wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    357                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    358                   hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    359               END DO 
    360             END DO 
    361          END DO 
    362  
    363          DO jl = 1, jpl 
     335                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) ) 
     336                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) ) 
     337                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) ) 
     338                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) ) 
     339                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) ) 
     340                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) ) 
     341               END DO 
     342            END DO 
     343 
    364344            DO jk = 1, nlay_i 
    365345               DO jj = 1, jpj 
    366346                  DO ji = 1, jpi 
    367                      rswitch          = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 
    368                      zei              = zs0e(ji,jj,jk,jl)       
    369                      e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 
    370                      ! Update fluxes 
    371                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    372                   END DO !ji 
    373                END DO ! jj 
    374             END DO ! jk 
    375          END DO ! jl 
    376  
    377          !--- Thickness correction in case too high (clem) -------------------------------------------------------- 
    378          CALL lim_var_glo2eqv 
     347                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 
     348                  END DO 
     349               END DO 
     350            END DO 
     351         END DO 
     352!!gm & cr  
     353 
     354         ! --- diags --- 
     355         DO jj = 1, jpj 
     356            DO ji = 1, jpi 
     357               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 
     358               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 
     359 
     360               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice 
     361               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice 
     362               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 
     363            END DO 
     364         END DO 
     365 
     366         ! zap small areas 
     367         CALL lim_var_zapsmall 
     368 
     369         !--- Thickness correction in case too high -------------------------------------------------------- 
    379370         DO jl = 1, jpl 
    380371            DO jj = 1, jpj 
     
    382373 
    383374                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
     375 
     376                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     377                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     378                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     379                      
    384380                     zvi  = v_i  (ji,jj,jl) 
    385381                     zvs  = v_s  (ji,jj,jl) 
     
    387383                     zes  = e_s  (ji,jj,1,jl) 
    388384                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) ) 
    389                      zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)    
    390                      !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)    
    391                       
    392                      rswitch = 1._wp 
    393                      IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 
    394                         & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                           
    395                         ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 
    396                         rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 
    397                         a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 
    398                      ELSE 
    399                         ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 
    400                         rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 
    401                         a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 
     385 
     386                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)   
     387 
     388                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
     389                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 
     390 
     391                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 
     392                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 
     393 
     394                        ! small correction due to *rswitch for a_i 
     395                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl) 
     396                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl) 
     397                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl) 
     398                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl) 
     399                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 
     400 
     401                        ! Update mass fluxes 
     402                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
     403                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
     404                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
     405                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
     406                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 
     407 
    402408                     ENDIF 
    403409 
    404                      ! small correction due to *rswitch for a_i 
    405                      v_i  (ji,jj,jl) = rswitch * v_i  (ji,jj,jl) 
    406                      v_s  (ji,jj,jl) = rswitch * v_s  (ji,jj,jl) 
    407                      smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 
    408                      e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 
    409                      e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 
    410  
    411                      ! Update mass fluxes 
    412                      wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
    413                      wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    414                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    415                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    416                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    417410                  ENDIF 
     411 
    418412               END DO 
    419413            END DO 
    420414         END DO 
    421415         ! ------------------------------------------------- 
    422  
    423          ! --- diags --- 
    424          DO jj = 1, jpj 
    425             DO ji = 1, jpi 
    426                diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 
    427                diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 
    428  
    429                diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 
    430                diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 
    431             END DO 
    432          END DO 
     416          
     417         !-------------------------------------- 
     418         ! Impose a_i < amax in mono-category 
     419         !-------------------------------------- 
     420         ! 
     421         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 
     422            DO jj = 1, jpj 
     423               DO ji = 1, jpi 
     424                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax ) 
     425               END DO 
     426            END DO 
     427         ENDIF 
    433428 
    434429         ! --- agglomerate variables ----------------- 
     
    436431         vt_s (:,:) = 0._wp 
    437432         at_i (:,:) = 0._wp 
    438          ! 
    439433         DO jl = 1, jpl 
    440434            DO jj = 1, jpj 
    441435               DO ji = 1, jpi 
    442                   ! 
    443                   vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 
    444                   vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 
    445                   at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    446                END DO 
    447             END DO 
    448          END DO 
    449          ! ------------------------------------------------- 
    450  
    451          ! open water 
     436                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) 
     437                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) 
     438                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 
     439               END DO 
     440            END DO 
     441         END DO 
     442 
     443         ! --- open water = 1 if at_i=0 -------------------------------- 
    452444         DO jj = 1, jpj 
    453445            DO ji = 1, jpi 
    454                ! open water = 1 if at_i=0 
    455446               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 
    456                ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj) 
     447               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 
    457448            END DO 
    458449         END DO       
     
    463454      ENDIF 
    464455 
    465       IF(ln_ctl) THEN   ! Control print 
    466          CALL prt_ctl_info(' ') 
    467          CALL prt_ctl_info(' - Cell values : ') 
    468          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    469          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :') 
    470          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :') 
    471          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :') 
    472          CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :') 
    473          DO jl = 1, jpl 
    474             CALL prt_ctl_info(' ') 
    475             CALL prt_ctl_info(' - Category : ', ivar1=jl) 
    476             CALL prt_ctl_info('   ~~~~~~~~~~') 
    477             CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ') 
    478             CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ') 
    479             CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ') 
    480             CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ') 
    481             CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ') 
    482             CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ') 
    483             CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ') 
    484             CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ') 
    485             CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ') 
    486             CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ') 
    487             DO jk = 1, nlay_i 
    488                CALL prt_ctl_info(' ') 
    489                CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    490                CALL prt_ctl_info('   ~~~~~~~') 
    491                CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ') 
    492                CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ') 
    493             END DO 
    494          END DO 
    495       ENDIF 
     456      ! ------------------------------------------------- 
     457      ! control prints 
     458      ! ------------------------------------------------- 
     459      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    496460      ! 
    497       CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 
    498       CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 
    499       CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 
    500  
    501       CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax )   ! clem 
     461      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold ) 
     462      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     463      CALL wrk_dealloc( jpi,jpj,1,          z0opw ) 
     464      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     465      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold ) 
    502466      ! 
    503467      IF( nn_timing == 1 )  CALL timing_stop('limtrp') 
     468 
    504469   END SUBROUTINE lim_trp 
    505470 
     
    512477   END SUBROUTINE lim_trp 
    513478#endif 
    514  
    515479   !!====================================================================== 
    516480END MODULE limtrp 
Note: See TracChangeset for help on using the changeset viewer.