Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90	(revision 11331)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90	(revision 11332)
@@ -773,37 +773,83 @@
    !! ** Purpose :  converting N-cat ice to jpl ice categories
    !!-------------------------------------------------------------------
-   SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting 1-cat ice to 1 ice category
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
-      REAL(wp), DIMENSION(:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
-      !!-------------------------------------------------------------------
-      zh_i(:) = zhti(:)
-      zh_s(:) = zhts(:)
-      za_i(:) = zati(:)
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
+      REAL(wp), DIMENSION(:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
+      REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      !!-------------------------------------------------------------------
+      ! == thickness and concentration == !
+      ph_i(:) = phti(:)
+      ph_s(:) = phts(:)
+      pa_i(:) = pati(:)
+      !
+      ! == temperature and salinity == !
+      IF( PRESENT( pt_i  ) )   pt_i (:) = ptmi (:)
+      IF( PRESENT( pt_s  ) )   pt_s (:) = ptms (:)
+      IF( PRESENT( pt_su ) )   pt_su(:) = ptmsu(:)
+      IF( PRESENT( ps_i  ) )   ps_i (:) = psmi (:)
+      
    END SUBROUTINE ice_var_itd_1c1c
 
-   SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting N-cat ice to 1 ice category
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
-      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
-      !!-------------------------------------------------------------------
-      !
-      za_i(:) = SUM( zati(:,:), dim=2 )
-      !
-      WHERE( za_i(:) /= 0._wp )
-         zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:)
-         zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:)
-      ELSEWHERE
-         zh_i(:) = 0._wp
-         zh_s(:) = 0._wp
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
+      REAL(wp), DIMENSION(:)  , INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
+      !
+      INTEGER ::   idim  
+      !!-------------------------------------------------------------------
+      !
+      idim = SIZE( phti, 1 )
+      !
+      ! == thickness and concentration == !
+      ALLOCATE( z1_ai(idim) )
+      !
+      pa_i(:) = SUM( pati(:,:), dim=2 )
+
+      WHERE( ( pa_i(:) ) /= 0._wp )   ;   z1_ai(:) = 1._wp / pa_i(:)
+      ELSEWHERE                       ;   z1_ai(:) = 0._wp
       END WHERE
+
+      ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
+      ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
+      !
+      ! == temperature and salinity == !
+      IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
+         !
+         ALLOCATE( z1_vi(idim), z1_vs(idim) )
+         !
+         WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
+         ELSEWHERE                                 ;   z1_vi(:) = 0._wp
+         END WHERE
+         WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
+         ELSEWHERE                                 ;   z1_vs(:) = 0._wp
+         END WHERE
+         !
+         IF( PRESENT( pt_i  ) )   pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+         IF( PRESENT( pt_s  ) )   pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
+         IF( PRESENT( pt_su ) )   pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
+         IF( PRESENT( ps_i  ) )   ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+         !
+         DEALLOCATE( z1_vi, z1_vs )
+         !
+      ENDIF
+      !
+      DEALLOCATE( z1_ai )
       !
    END SUBROUTINE ice_var_itd_Nc1c
    
-   SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
       !!-------------------------------------------------------------------
       !!
@@ -826,7 +872,7 @@
       !!               4) Iterate until ok (SUM(itest(:) = 4)
       !!
-      !! ** Arguments : zhti: 1-cat ice thickness
-      !!                zhts: 1-cat snow depth
-      !!                zati: 1-cat ice concentration
+      !! ** Arguments : phti: 1-cat ice thickness
+      !!                phts: 1-cat snow depth
+      !!                pati: 1-cat ice concentration
       !!
       !! ** Output    : jpl-cat 
@@ -834,36 +880,34 @@
       !!  (Example of application: BDY forcings when input are cell averaged)  
       !!-------------------------------------------------------------------
-      INTEGER  :: ji, jk, jl             ! dummy loop indices
-      INTEGER  :: idim, i_fill, jl0  
-      REAL(wp) :: zarg, zV, zconv, zdh, zdv
-      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
-      INTEGER , DIMENSION(4)                  ::   itest
-      !!-------------------------------------------------------------------
-      !
-      ! ----------------------------------------
+      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
+      REAL(wp), DIMENSION(:)  , INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
+      REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      !
+      INTEGER , DIMENSION(4) ::   itest
+      INTEGER  ::   ji, jk, jl
+      INTEGER  ::   idim, i_fill, jl0  
+      REAL(wp) ::   zarg, zV, zconv, zdh, zdv
+      !!-------------------------------------------------------------------
+      !
+      ! == thickness and concentration == !
       ! distribution over the jpl ice categories
-      ! ----------------------------------------
-      ! a gaussian distribution for ice concentration is used
-      ! then we check whether the distribution fullfills
-      ! volume and area conservation, positivity and ice categories bounds
-      idim = SIZE( zhti , 1 )
-      zh_i(1:idim,1:jpl) = 0._wp
-      zh_s(1:idim,1:jpl) = 0._wp
-      za_i(1:idim,1:jpl) = 0._wp
-      !
-      IF( jpl == 1 ) THEN
-         CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
-         RETURN
-      ENDIF
+      !    a gaussian distribution for ice concentration is used
+      !    then we check whether the distribution fullfills
+      !    volume and area conservation, positivity and ice categories bounds
+      idim = SIZE( phti , 1 )
+      !
+      ph_i(1:idim,1:jpl) = 0._wp
+      ph_s(1:idim,1:jpl) = 0._wp
+      pa_i(1:idim,1:jpl) = 0._wp
       !
       DO ji = 1, idim
          !
-         IF( zhti(ji) > 0._wp ) THEN
+         IF( phti(ji) > 0._wp ) THEN
             !
             ! find which category (jl0) the input ice thickness falls into
             jl0 = jpl
             DO jl = 1, jpl
-               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+               IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN
                   jl0 = jl
                   CYCLE
@@ -877,30 +921,30 @@
                i_fill = i_fill - 1
                !
-               zh_i(ji,1:jpl) = 0._wp
-               za_i(ji,1:jpl) = 0._wp
+               ph_i(ji,1:jpl) = 0._wp
+               pa_i(ji,1:jpl) = 0._wp
                itest(:)       = 0      
                !
                IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
-                  zh_i(ji,1) = zhti(ji)
-                  za_i (ji,1) = zati (ji)
+                  ph_i(ji,1) = phti(ji)
+                  pa_i (ji,1) = pati (ji)
                ELSE                         !-- case ice is thicker: fill categories >1
                   ! thickness
                   DO jl = 1, i_fill - 1
-                     zh_i(ji,jl) = hi_mean(jl)
+                     ph_i(ji,jl) = hi_mean(jl)
                   END DO
                   !
                   ! concentration
-                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
+                  pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl))
                   DO jl = 1, i_fill - 1
                      IF ( jl /= jl0 ) THEN
-                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
-                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
+                        zarg        = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp )
+                        pa_i(ji,jl) =   pa_i (ji,jl0) * EXP(-zarg**2)
                      ENDIF
                   END DO
                   !
                   ! last category
-                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
-                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
-                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
+                  pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) )
+                  zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) )
+                  ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 
                   !
                   ! correction if concentration of upper cat is greater than lower cat
@@ -908,9 +952,9 @@
                   IF ( jl0 /= jpl ) THEN
                      DO jl = jpl, jl0+1, -1
-                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
-                           zdv = zh_i(ji,jl) * za_i(ji,jl)
-                           zh_i(ji,jl    ) = 0._wp
-                           za_i (ji,jl    ) = 0._wp
-                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
+                        IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN
+                           zdv = ph_i(ji,jl) * pa_i(ji,jl)
+                           ph_i(ji,jl    ) = 0._wp
+                           pa_i (ji,jl    ) = 0._wp
+                           pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 )
                         END IF
                      END DO
@@ -920,15 +964,15 @@
                !
                ! Compatibility tests
-               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
+               zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 
                IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
                !
-               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
+               zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) )
                IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
                !
-               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
+               IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
                !
                itest(4) = 1
                DO jl = 1, i_fill
-                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
+                  IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
                END DO
                !                                         !----------------------------
@@ -938,22 +982,45 @@
       END DO
 
-      ! Add Snow in each category where za_i is not 0
+      ! Add Snow in each category where pa_i is not 0
       DO jl = 1, jpl
          DO ji = 1, idim
-            IF( za_i(ji,jl) > 0._wp ) THEN
-               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
+            IF( pa_i(ji,jl) > 0._wp ) THEN
+               ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) )
                ! In case snow load is in excess that would lead to transformation from snow to ice
                ! Then, transfer the snow excess into the ice (different from icethd_dh)
-               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
+               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 
                ! recompute h_i, h_s avoiding out of bounds values
-               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
-               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
+               ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
+               ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
             ENDIF
          END DO
       END DO
       !
+      ! == temperature and salinity == !
+      IF( PRESENT( pt_i  ) ) THEN
+         DO jl = 1, jpl
+            pt_i(:,jl) = ptmi (:)
+         END DO
+      ENDIF
+      IF( PRESENT( pt_s  ) ) THEN
+         DO jl = 1, jpl
+            pt_s (:,jl) = ptms (:)
+         END DO
+      ENDIF
+      IF( PRESENT( pt_su ) ) THEN
+         DO jl = 1, jpl
+            pt_su(:,jl) = ptmsu(:)
+         END DO
+      ENDIF
+      IF( PRESENT( ps_i  ) ) THEN
+         DO jl = 1, jpl
+            ps_i (:,jl) = psmi (:)
+         END DO
+      ENDIF
+      !
    END SUBROUTINE ice_var_itd_1cMc
 
-   SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,       ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
       !!-------------------------------------------------------------------
       !!
@@ -976,7 +1043,7 @@
       !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
       !!
-      !! ** Arguments : zhti: N-cat ice thickness
-      !!                zhts: N-cat snow depth
-      !!                zati: N-cat ice concentration
+      !! ** Arguments : phti: N-cat ice thickness
+      !!                phts: N-cat snow depth
+      !!                pati: N-cat ice concentration
       !!
       !! ** Output    : jpl-cat 
@@ -984,33 +1051,44 @@
       !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)  
       !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
+      REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal
+      REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal
+      !
+      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
+      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin
+      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp
+      !
+      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
+      INTEGER  ::   ji, jl, jl1, jl2
       INTEGER  ::   idim, icat  
-      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
-      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
-      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
-      !!-------------------------------------------------------------------
-      !
-      idim = SIZE( zhti, 1 )
-      icat = SIZE( zhti, 2 )
+      !!-------------------------------------------------------------------
+      !
+      idim = SIZE( phti, 1 )
+      icat = SIZE( phti, 2 )
+      !
+      ! == thickness and concentration == !
       !                                 ! ---------------------- !
       IF( icat == jpl ) THEN            ! input cat = output cat !
          !                              ! ---------------------- !
-         zh_i(:,:) = zhti(:,:)
-         zh_s(:,:) = zhts(:,:)
-         za_i(:,:) = zati(:,:)
+         ph_i(:,:) = phti(:,:)
+         ph_s(:,:) = phts(:,:)
+         pa_i(:,:) = pati(:,:)
+         !
+         ! == temperature and salinity == !
+         IF( PRESENT( pt_i  ) )   pt_i (:,:) = ptmi (:,:)
+         IF( PRESENT( pt_s  ) )   pt_s (:,:) = ptms (:,:)
+         IF( PRESENT( pt_su ) )   pt_su(:,:) = ptmsu(:,:)
+         IF( PRESENT( ps_i  ) )   ps_i (:,:) = psmi (:,:)
          !                              ! ---------------------- !
-      ELSEIF( icat == 1 ) THEN          ! specific case if N = 1 !
+      ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
          !                              ! ---------------------- !
-         !
-         CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i )
-         !
+         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), &
+            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) )
          !                              ! ---------------------- !
-      ELSEIF( jpl == 1 ) THEN           ! specific case if M = 1 !
+      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
          !                              ! ---------------------- !
-         !
-         CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) )
-         !
+         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), &
+            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )         
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !
@@ -1021,7 +1099,7 @@
 
          ! --- initialize output fields to 0 --- !
-         zh_i(1:idim,1:jpl) = 0._wp
-         zh_s(1:idim,1:jpl) = 0._wp
-         za_i(1:idim,1:jpl) = 0._wp
+         ph_i(1:idim,1:jpl) = 0._wp
+         ph_s(1:idim,1:jpl) = 0._wp
+         pa_i(1:idim,1:jpl) = 0._wp
          !
          ! --- fill the categories --- !
@@ -1033,9 +1111,9 @@
             DO jl2 = 1, icat
                DO ji = 1, idim
-                  IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
+                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
                      ! fill the right category
-                     zh_i(ji,jl1) = zhti(ji,jl2)
-                     zh_s(ji,jl1) = zhts(ji,jl2)
-                     za_i(ji,jl1) = zati(ji,jl2)
+                     ph_i(ji,jl1) = phti(ji,jl2)
+                     ph_s(ji,jl1) = phts(ji,jl2)
+                     pa_i(ji,jl1) = pati(ji,jl2)
                      ! record categories that are filled
                      jlmax(ji) = MAX( jlmax(ji), jl1 )
@@ -1054,15 +1132,15 @@
             IF( jl1 > 1 ) THEN
                ! fill the lower cat (jl1-1)
-               za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
-               zh_i(ji,jl1-1) = hi_mean(jl1-1)
+               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
+               ph_i(ji,jl1-1) = hi_mean(jl1-1)
                ! remove from cat jl1
-               za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
+               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
             ENDIF
             IF( jl2 < jpl ) THEN
                ! fill the upper cat (jl2+1)
-               za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
-               zh_i(ji,jl2+1) = hi_mean(jl2+1)
+               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
+               ph_i(ji,jl2+1) = hi_mean(jl2+1)
                ! remove from cat jl2
-               za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
+               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
             ENDIF
          END DO
@@ -1074,9 +1152,9 @@
                IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
                   ! fill high
-                  za_i(ji,jl) = ztrans * za_i(ji,jl-1)
-                  zh_i(ji,jl) = hi_mean(jl)
+                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
+                  ph_i(ji,jl) = hi_mean(jl)
                   jlfil(ji,jl) = jl
                   ! remove low
-                  za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
+                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
                ENDIF
             END DO
@@ -1088,9 +1166,9 @@
                IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
                   ! fill low
-                  za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
-                  zh_i(ji,jl) = hi_mean(jl) 
+                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
+                  ph_i(ji,jl) = hi_mean(jl) 
                   jlfil2(ji,jl) = jl
                   ! remove high
-                  za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
+                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
                ENDIF
             END DO
@@ -1099,4 +1177,50 @@
          DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
          DEALLOCATE( jlmin, jlmax )
+         !
+         ! == temperature and salinity == !
+         !
+         IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
+            !
+            ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
+            !
+            WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
+            ELSEWHERE                                               ;   z1_ai(:) = 0._wp
+            END WHERE
+            WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
+            ELSEWHERE                                               ;   z1_vi(:) = 0._wp
+            END WHERE
+            WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
+            ELSEWHERE                                               ;   z1_vs(:) = 0._wp
+            END WHERE
+            !
+            ! fill all the categories with the same value
+            IF( PRESENT( pt_i  ) ) THEN
+               ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+               DO jl = 1, jpl
+                  pt_i (:,jl) = ztmp(:)
+               END DO
+            ENDIF
+            IF( PRESENT( pt_s  ) ) THEN
+               ztmp(:) =  SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
+               DO jl = 1, jpl
+                  pt_s (:,jl) = ztmp(:)
+               END DO
+            ENDIF
+            IF( PRESENT( pt_su ) ) THEN
+               ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
+               DO jl = 1, jpl
+                  pt_su(:,jl) = ztmp(:)
+               END DO
+            ENDIF
+            IF( PRESENT( ps_i  ) ) THEN
+               ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
+               DO jl = 1, jpl
+                  ps_i (:,jl) = ztmp(:)
+               END DO
+            ENDIF
+            !
+            DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
+            !
+         ENDIF
          !
       ENDIF
