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 12178 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icevar.F90 – NEMO

Ignore:
Timestamp:
2019-12-11T12:02:38+01:00 (4 years ago)
Author:
agn
Message:

updated trunk to v 11653

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icevar.F90

    r10994 r12178  
    3232   !!                        - vt_s(jpi,jpj) 
    3333   !!                        - at_i(jpi,jpj) 
     34   !!                        - st_i(jpi,jpj) 
    3435   !!                        - et_s(jpi,jpj)  total snow heat content 
    3536   !!                        - et_i(jpi,jpj)  total ice thermal content  
     
    4647   !!   ice_var_zapneg    : remove negative ice fields 
    4748   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    48    !!   ice_var_itd       : convert 1-cat to jpl-cat 
    49    !!   ice_var_itd2      : convert N-cat to jpl-cat 
    5049   !!   ice_var_bv        : brine volume 
    5150   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature 
    5251   !!   ice_var_sshdyn    : compute equivalent ssh in lead 
     52   !!   ice_var_itd       : convert N-cat to M-cat 
    5353   !!---------------------------------------------------------------------- 
    5454   USE dom_oce        ! ocean space and time domain 
     
    7373   PUBLIC   ice_var_zapneg 
    7474   PUBLIC   ice_var_roundoff 
    75    PUBLIC   ice_var_itd 
    76    PUBLIC   ice_var_itd2 
    7775   PUBLIC   ice_var_bv            
    7876   PUBLIC   ice_var_enthalpy            
    7977   PUBLIC   ice_var_sshdyn 
     78   PUBLIC   ice_var_itd 
     79 
     80   INTERFACE ice_var_itd 
     81      MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc 
     82   END INTERFACE 
    8083 
    8184   !!---------------------------------------------------------------------- 
     
    101104      ! 
    102105      !                                      ! integrated values 
    103       vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 ) 
    104       vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 ) 
    105       at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 ) 
    106       et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
    107       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
     106      vt_i(:,:) =       SUM( v_i (:,:,:)           , dim=3 ) 
     107      vt_s(:,:) =       SUM( v_s (:,:,:)           , dim=3 ) 
     108      st_i(:,:) =       SUM( sv_i(:,:,:)           , dim=3 ) 
     109      at_i(:,:) =       SUM( a_i (:,:,:)           , dim=3 ) 
     110      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 
     111      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 
    108112      ! 
    109113      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
     
    135139         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    136140         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
    137          sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:) 
     141         sm_i (:,:) =      st_i(:,:)                          * z1_vt_i(:,:) 
    138142         ! 
    139143         tm_i(:,:) = 0._wp 
     
    155159            tm_s (:,:) = rt0 
    156160         END WHERE 
    157  
     161         ! 
     162         !                           ! mean melt pond depth 
     163         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 
     164         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp 
     165         END WHERE          
     166         ! 
    158167         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 
     168         ! 
    159169      ENDIF 
    160170      ! 
     
    260270      ! 
    261271      ! integrated values  
    262       vt_i (:,:) = SUM( v_i, dim=3 ) 
    263       vt_s (:,:) = SUM( v_s, dim=3 ) 
    264       at_i (:,:) = SUM( a_i, dim=3 ) 
     272      vt_i (:,:) = SUM( v_i , dim=3 ) 
     273      vt_s (:,:) = SUM( v_s , dim=3 ) 
     274      at_i (:,:) = SUM( a_i , dim=3 ) 
    265275      ! 
    266276   END SUBROUTINE ice_var_glo2eqv 
     
    530540 
    531541      ! to be sure that at_i is the sum of a_i(jl) 
    532       at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 
    533       vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 
     542      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) 
     543      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) 
     544!!clem add? 
     545!      vt_s (:,:) = SUM( v_s (:,:,:), dim=3 ) 
     546!      st_i (:,:) = SUM( sv_i(:,:,:), dim=3 ) 
     547!      et_s(:,:)  = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 
     548!      et_i(:,:)  = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 
     549!!clem 
    534550 
    535551      ! open water = 1 if at_i=0 
     
    649665      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    650666      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    651       IF ( ln_pnd_H12 ) THEN 
     667      IF( ln_pnd_H12 ) THEN 
    652668         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    653669         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     
    656672   END SUBROUTINE ice_var_roundoff 
    657673    
    658     
    659    SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    660       !!------------------------------------------------------------------- 
    661       !!                ***  ROUTINE ice_var_itd   *** 
    662       !! 
    663       !! ** Purpose :  converting 1-cat ice to multiple ice categories 
    664       !! 
    665       !!                  ice thickness distribution follows a gaussian law 
    666       !!               around the concentration of the most likely ice thickness 
    667       !!                           (similar as iceistate.F90) 
    668       !! 
    669       !! ** Method:   Iterative procedure 
    670       !!                 
    671       !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
    672       !! 
    673       !!               2) Check whether the distribution conserves area and volume, positivity and 
    674       !!                  category boundaries 
    675       !!               
    676       !!               3) If not (input ice is too thin), the last category is empty and 
    677       !!                  the number of categories is reduced (jpl-1) 
    678       !! 
    679       !!               4) Iterate until ok (SUM(itest(:) = 4) 
    680       !! 
    681       !! ** Arguments : zhti: 1-cat ice thickness 
    682       !!                zhts: 1-cat snow depth 
    683       !!                zati: 1-cat ice concentration 
    684       !! 
    685       !! ** Output    : jpl-cat  
    686       !! 
    687       !!  (Example of application: BDY forcings when input are cell averaged)   
    688       !!------------------------------------------------------------------- 
    689       INTEGER  :: ji, jk, jl             ! dummy loop indices 
    690       INTEGER  :: idim, i_fill, jl0   
    691       REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    692       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables 
    693       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    694       INTEGER , DIMENSION(4)                  ::   itest 
    695       !!------------------------------------------------------------------- 
    696       ! 
    697       ! ---------------------------------------- 
    698       ! distribution over the jpl ice categories 
    699       ! ---------------------------------------- 
    700       ! a gaussian distribution for ice concentration is used 
    701       ! then we check whether the distribution fullfills 
    702       ! volume and area conservation, positivity and ice categories bounds 
    703       idim = SIZE( zhti , 1 ) 
    704       zh_i(1:idim,1:jpl) = 0._wp 
    705       zh_s(1:idim,1:jpl) = 0._wp 
    706       za_i(1:idim,1:jpl) = 0._wp 
    707       ! 
    708       DO ji = 1, idim 
    709          ! 
    710          IF( zhti(ji) > 0._wp ) THEN 
    711             ! 
    712             ! find which category (jl0) the input ice thickness falls into 
    713             jl0 = jpl 
    714             DO jl = 1, jpl 
    715                IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
    716                   jl0 = jl 
    717                   CYCLE 
    718                ENDIF 
    719             END DO 
    720             ! 
    721             itest(:) = 0 
    722             i_fill   = jpl + 1                                            !------------------------------------ 
    723             DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    724                !                                                          !------------------------------------ 
    725                i_fill = i_fill - 1 
    726                ! 
    727                zh_i(ji,1:jpl) = 0._wp 
    728                za_i(ji,1:jpl) = 0._wp 
    729                itest(:)       = 0       
    730                ! 
    731                IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    732                   zh_i(ji,1) = zhti(ji) 
    733                   za_i (ji,1) = zati (ji) 
    734                ELSE                         !-- case ice is thicker: fill categories >1 
    735                   ! thickness 
    736                   DO jl = 1, i_fill - 1 
    737                      zh_i(ji,jl) = hi_mean(jl) 
    738                   END DO 
    739                   ! 
    740                   ! concentration 
    741                   za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 
    742                   DO jl = 1, i_fill - 1 
    743                      IF ( jl /= jl0 ) THEN 
    744                         zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    745                         za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
    746                      ENDIF 
    747                   END DO 
    748                   ! 
    749                   ! last category 
    750                   za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) ) 
    751                   zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) ) 
    752                   zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )  
    753                   ! 
    754                   ! correction if concentration of upper cat is greater than lower cat 
    755                   !    (it should be a gaussian around jl0 but sometimes it is not) 
    756                   IF ( jl0 /= jpl ) THEN 
    757                      DO jl = jpl, jl0+1, -1 
    758                         IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 
    759                            zdv = zh_i(ji,jl) * za_i(ji,jl) 
    760                            zh_i(ji,jl    ) = 0._wp 
    761                            za_i (ji,jl    ) = 0._wp 
    762                            za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 
    763                         END IF 
    764                      END DO 
    765                   ENDIF 
    766                   ! 
    767                ENDIF 
    768                ! 
    769                ! Compatibility tests 
    770                zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )  
    771                IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation 
    772                ! 
    773                zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 
    774                IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation 
    775                ! 
    776                IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ? 
    777                ! 
    778                itest(4) = 1 
    779                DO jl = 1, i_fill 
    780                   IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations 
    781                END DO 
    782                !                                         !---------------------------- 
    783             END DO                                       ! end iteration on categories 
    784             !                                            !---------------------------- 
    785          ENDIF 
    786       END DO 
    787  
    788       ! Add Snow in each category where za_i is not 0 
    789       DO jl = 1, jpl 
    790          DO ji = 1, idim 
    791             IF( za_i(ji,jl) > 0._wp ) THEN 
    792                zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
    793                ! In case snow load is in excess that would lead to transformation from snow to ice 
    794                ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    795                zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )  
    796                ! recompute h_i, h_s avoiding out of bounds values 
    797                zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh ) 
    798                zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos ) 
    799             ENDIF 
    800          END DO 
    801       END DO 
    802       ! 
    803    END SUBROUTINE ice_var_itd 
    804  
    805  
    806    SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    807       !!------------------------------------------------------------------- 
    808       !!                ***  ROUTINE ice_var_itd2   *** 
    809       !! 
    810       !! ** Purpose :  converting N-cat ice to jpl ice categories 
    811       !! 
    812       !!                  ice thickness distribution follows a gaussian law 
    813       !!               around the concentration of the most likely ice thickness 
    814       !!                           (similar as iceistate.F90) 
    815       !! 
    816       !! ** Method:   Iterative procedure 
    817       !!                 
    818       !!               1) Fill ice cat that correspond to input thicknesses 
    819       !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled 
    820       !! 
    821       !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    822        !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    823       !!               
    824       !!               3) Expand the filling to the empty cat between jlmin and jlmax  
    825       !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 
    826       !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 
    827       !! 
    828       !! ** Arguments : zhti: N-cat ice thickness 
    829       !!                zhts: N-cat snow depth 
    830       !!                zati: N-cat ice concentration 
    831       !! 
    832       !! ** Output    : jpl-cat  
    833       !! 
    834       !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)   
    835       !!------------------------------------------------------------------- 
    836       INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices 
    837       INTEGER  ::   idim, icat   
    838       REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
    839       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    840       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    841       INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2 
    842       INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin 
    843       !!------------------------------------------------------------------- 
    844       ! 
    845       idim = SIZE( zhti, 1 ) 
    846       icat = SIZE( zhti, 2 ) 
    847       ! 
    848       ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays 
    849       ALLOCATE( jlmin(idim), jlmax(idim) ) 
    850  
    851       ! --- initialize output fields to 0 --- ! 
    852       zh_i(1:idim,1:jpl) = 0._wp 
    853       zh_s(1:idim,1:jpl) = 0._wp 
    854       za_i(1:idim,1:jpl) = 0._wp 
    855       ! 
    856       ! --- fill the categories --- ! 
    857       !     find where cat-input = cat-output and fill cat-output fields   
    858       jlmax(:) = 0 
    859       jlmin(:) = 999 
    860       jlfil(:,:) = 0 
    861       DO jl1 = 1, jpl 
    862          DO jl2 = 1, icat 
    863             DO ji = 1, idim 
    864                IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
    865                   ! fill the right category 
    866                   zh_i(ji,jl1) = zhti(ji,jl2) 
    867                   zh_s(ji,jl1) = zhts(ji,jl2) 
    868                   za_i(ji,jl1) = zati(ji,jl2) 
    869                   ! record categories that are filled 
    870                   jlmax(ji) = MAX( jlmax(ji), jl1 ) 
    871                   jlmin(ji) = MIN( jlmin(ji), jl1 ) 
    872                   jlfil(ji,jl1) = jl1 
    873                ENDIF 
    874             END DO 
    875          END DO 
    876       END DO 
    877       ! 
    878       ! --- fill the gaps between categories --- !   
    879       !     transfer from categories filled at the previous step to the empty ones in between 
    880       DO ji = 1, idim 
    881          jl1 = jlmin(ji) 
    882          jl2 = jlmax(ji) 
    883          IF( jl1 > 1 ) THEN 
    884             ! fill the lower cat (jl1-1) 
    885             za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 
    886             zh_i(ji,jl1-1) = hi_mean(jl1-1) 
    887             ! remove from cat jl1 
    888             za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
    889          ENDIF 
    890          IF( jl2 < jpl ) THEN 
    891             ! fill the upper cat (jl2+1) 
    892             za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 
    893             zh_i(ji,jl2+1) = hi_mean(jl2+1) 
    894             ! remove from cat jl2 
    895             za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
    896          ENDIF 
    897       END DO 
    898       ! 
    899       jlfil2(:,:) = jlfil(:,:)  
    900       ! fill categories from low to high 
    901       DO jl = 2, jpl-1 
    902          DO ji = 1, idim 
    903             IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
    904                ! fill high 
    905                za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
    906                zh_i(ji,jl) = hi_mean(jl) 
    907                jlfil(ji,jl) = jl 
    908                ! remove low 
    909                za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
    910             ENDIF 
    911          END DO 
    912       END DO 
    913       ! 
    914       ! fill categories from high to low 
    915       DO jl = jpl-1, 2, -1 
    916          DO ji = 1, idim 
    917             IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
    918                ! fill low 
    919                za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 
    920                zh_i(ji,jl) = hi_mean(jl)  
    921                jlfil2(ji,jl) = jl 
    922                ! remove high 
    923                za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
    924             ENDIF 
    925          END DO 
    926       END DO 
    927       ! 
    928       DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
    929       DEALLOCATE( jlmin, jlmax ) 
    930       ! 
    931    END SUBROUTINE ice_var_itd2 
    932  
    933674 
    934675   SUBROUTINE ice_var_bv 
     
    992733   END SUBROUTINE ice_var_enthalpy 
    993734 
     735    
    994736   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 
    995737      !!--------------------------------------------------------------------- 
     
    1038780   END FUNCTION ice_var_sshdyn 
    1039781 
     782    
     783   !!------------------------------------------------------------------- 
     784   !!                ***  INTERFACE ice_var_itd   *** 
     785   !! 
     786   !! ** Purpose :  converting N-cat ice to jpl ice categories 
     787   !!------------------------------------------------------------------- 
     788   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     789      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     790      !!------------------------------------------------------------------- 
     791      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     792      !!------------------------------------------------------------------- 
     793      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     794      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     795      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     796      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     797      !!------------------------------------------------------------------- 
     798      ! == thickness and concentration == ! 
     799      ph_i(:) = phti(:) 
     800      ph_s(:) = phts(:) 
     801      pa_i(:) = pati(:) 
     802      ! 
     803      ! == temperature and salinity and ponds == ! 
     804      pt_i (:) = ptmi (:) 
     805      pt_s (:) = ptms (:) 
     806      pt_su(:) = ptmsu(:) 
     807      ps_i (:) = psmi (:) 
     808      pa_ip(:) = patip(:) 
     809      ph_ip(:) = phtip(:) 
     810       
     811   END SUBROUTINE ice_var_itd_1c1c 
     812 
     813   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     814      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     815      !!------------------------------------------------------------------- 
     816      !! ** Purpose :  converting N-cat ice to 1 ice category 
     817      !!------------------------------------------------------------------- 
     818      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     819      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     820      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     821      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     822      ! 
     823      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
     824      ! 
     825      INTEGER ::   idim   
     826      !!------------------------------------------------------------------- 
     827      ! 
     828      idim = SIZE( phti, 1 ) 
     829      ! 
     830      ! == thickness and concentration == ! 
     831      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 
     832      ! 
     833      pa_i(:) = SUM( pati(:,:), dim=2 ) 
     834 
     835      WHERE( ( pa_i(:) ) /= 0._wp )   ;   z1_ai(:) = 1._wp / pa_i(:) 
     836      ELSEWHERE                       ;   z1_ai(:) = 0._wp 
     837      END WHERE 
     838 
     839      ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 
     840      ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 
     841      ! 
     842      ! == temperature and salinity == ! 
     843      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 
     844      ELSEWHERE                                 ;   z1_vi(:) = 0._wp 
     845      END WHERE 
     846      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 
     847      ELSEWHERE                                 ;   z1_vs(:) = 0._wp 
     848      END WHERE 
     849      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     850      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     851      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     852      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     853 
     854      ! == ponds == ! 
     855      pa_ip(:) = SUM( patip(:,:), dim=2 ) 
     856      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     857      ELSEWHERE                    ;   ph_ip(:) = 0._wp 
     858      END WHERE 
     859      ! 
     860      DEALLOCATE( z1_ai, z1_vi, z1_vs ) 
     861      ! 
     862   END SUBROUTINE ice_var_itd_Nc1c 
     863    
     864   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     865      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     866      !!------------------------------------------------------------------- 
     867      !! 
     868      !! ** Purpose :  converting 1-cat ice to jpl ice categories 
     869      !! 
     870      !! 
     871      !! ** Method:   ice thickness distribution follows a gamma function from Abraham et al. (2015) 
     872      !!              it has the property of conserving total concentration and volume 
     873      !!               
     874      !! 
     875      !! ** Arguments : phti: 1-cat ice thickness 
     876      !!                phts: 1-cat snow depth 
     877      !!                pati: 1-cat ice concentration 
     878      !! 
     879      !! ** Output    : jpl-cat  
     880      !! 
     881      !!  Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 
     882      !!               Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 
     883      !!               Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614  
     884      !!------------------------------------------------------------------- 
     885      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     886      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     887      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     888      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     889      ! 
     890      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti 
     891      INTEGER  ::   ji, jk, jl 
     892      INTEGER  ::   idim 
     893      REAL(wp) ::   zv, zdh 
     894      !!------------------------------------------------------------------- 
     895      ! 
     896      idim = SIZE( phti , 1 ) 
     897      ! 
     898      ph_i(1:idim,1:jpl) = 0._wp 
     899      ph_s(1:idim,1:jpl) = 0._wp 
     900      pa_i(1:idim,1:jpl) = 0._wp 
     901      ! 
     902      ALLOCATE( z1_hti(idim) ) 
     903      WHERE( phti(:) /= 0._wp )   ;   z1_hti(:) = 1._wp / phti(:) 
     904      ELSEWHERE                   ;   z1_hti(:) = 0._wp 
     905      END WHERE 
     906      ! 
     907      ! == thickness and concentration == ! 
     908      ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 
     909      DO jl = 1, jpl-1 
     910         DO ji = 1, idim 
     911            ! 
     912            IF( phti(ji) > 0._wp ) THEN 
     913               ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 
     914               pa_i(ji,jl) = pati(ji) * z1_hti(ji) * (  ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 
     915                  &                                   - ( phti(ji) + 2.*hi_max(jl  ) ) * EXP( -2.*hi_max(jl  )*z1_hti(ji) ) ) 
     916               ! 
     917               ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 
     918               zv = pati(ji) * z1_hti(ji) * (  ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 
     919                  &                            * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 
     920                  &                          - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 
     921                  &                            * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 
     922               ! thickness 
     923               IF( pa_i(ji,jl) > epsi06 ) THEN 
     924                  ph_i(ji,jl) = zv / pa_i(ji,jl) 
     925               ELSE 
     926                  ph_i(ji,jl) = 0. 
     927                  pa_i(ji,jl) = 0. 
     928               ENDIF 
     929            ENDIF 
     930            ! 
     931         ENDDO 
     932      ENDDO 
     933      ! 
     934      ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 
     935      DO ji = 1, idim 
     936         ! 
     937         IF( phti(ji) > 0._wp ) THEN 
     938            ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 
     939            pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 
     940 
     941            ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 
     942            zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 
     943               &                         * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 
     944            ! thickness 
     945            IF( pa_i(ji,jpl) > epsi06 ) THEN 
     946               ph_i(ji,jpl) = zv / pa_i(ji,jpl) 
     947            else 
     948               ph_i(ji,jpl) = 0. 
     949               pa_i(ji,jpl) = 0. 
     950            ENDIF 
     951         ENDIF 
     952         ! 
     953      ENDDO 
     954      ! 
     955      ! Add Snow in each category where pa_i is not 0 
     956      DO jl = 1, jpl 
     957         DO ji = 1, idim 
     958            IF( pa_i(ji,jl) > 0._wp ) THEN 
     959               ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 
     960               ! In case snow load is in excess that would lead to transformation from snow to ice 
     961               ! Then, transfer the snow excess into the ice (different from icethd_dh) 
     962               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )  
     963               ! recompute h_i, h_s avoiding out of bounds values 
     964               ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 
     965               ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 
     966            ENDIF 
     967         END DO 
     968      END DO 
     969      ! 
     970      DEALLOCATE( z1_hti ) 
     971      ! 
     972      ! == temperature and salinity == ! 
     973      DO jl = 1, jpl 
     974         pt_i (:,jl) = ptmi (:) 
     975         pt_s (:,jl) = ptms (:) 
     976         pt_su(:,jl) = ptmsu(:) 
     977         ps_i (:,jl) = psmi (:) 
     978         ps_i (:,jl) = psmi (:)          
     979      END DO 
     980      ! 
     981      ! == ponds == ! 
     982      ALLOCATE( zfra(idim) ) 
     983      ! keep the same pond fraction atip/ati for each category 
     984      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:) 
     985      ELSEWHERE                   ;   zfra(:) = 0._wp 
     986      END WHERE 
     987      DO jl = 1, jpl 
     988         pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     989      END DO 
     990      ! keep the same v_ip/v_i ratio for each category 
     991      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     992      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     993      END WHERE 
     994      DO jl = 1, jpl 
     995         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     996         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     997         END WHERE 
     998      END DO 
     999      DEALLOCATE( zfra ) 
     1000      ! 
     1001   END SUBROUTINE ice_var_itd_1cMc 
     1002 
     1003   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     1004      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     1005      !!------------------------------------------------------------------- 
     1006      !! 
     1007      !! ** Purpose :  converting N-cat ice to jpl ice categories 
     1008      !! 
     1009      !!                  ice thickness distribution follows a gaussian law 
     1010      !!               around the concentration of the most likely ice thickness 
     1011      !!                           (similar as iceistate.F90) 
     1012      !! 
     1013      !! ** Method:   Iterative procedure 
     1014      !!                 
     1015      !!               1) Fill ice cat that correspond to input thicknesses 
     1016      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled 
     1017      !! 
     1018      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
     1019       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     1020      !!               
     1021      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
     1022      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 
     1023      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 
     1024      !! 
     1025      !! ** Arguments : phti: N-cat ice thickness 
     1026      !!                phts: N-cat snow depth 
     1027      !!                pati: N-cat ice concentration 
     1028      !! 
     1029      !! ** Output    : jpl-cat  
     1030      !! 
     1031      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)   
     1032      !!------------------------------------------------------------------- 
     1033      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     1034      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     1035      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     1036      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     1037      ! 
     1038      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
     1039      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin 
     1040      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra 
     1041      ! 
     1042      REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
     1043      INTEGER  ::   ji, jl, jl1, jl2 
     1044      INTEGER  ::   idim, icat   
     1045      !!------------------------------------------------------------------- 
     1046      ! 
     1047      idim = SIZE( phti, 1 ) 
     1048      icat = SIZE( phti, 2 ) 
     1049      ! 
     1050      ! == thickness and concentration == ! 
     1051      !                                 ! ---------------------- ! 
     1052      IF( icat == jpl ) THEN            ! input cat = output cat ! 
     1053         !                              ! ---------------------- ! 
     1054         ph_i(:,:) = phti(:,:) 
     1055         ph_s(:,:) = phts(:,:) 
     1056         pa_i(:,:) = pati(:,:) 
     1057         ! 
     1058         ! == temperature and salinity and ponds == ! 
     1059         pt_i (:,:) = ptmi (:,:) 
     1060         pt_s (:,:) = ptms (:,:) 
     1061         pt_su(:,:) = ptmsu(:,:) 
     1062         ps_i (:,:) = psmi (:,:) 
     1063         pa_ip(:,:) = patip(:,:) 
     1064         ph_ip(:,:) = phtip(:,:) 
     1065         !                              ! ---------------------- ! 
     1066      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
     1067         !                              ! ---------------------- ! 
     1068         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
     1069            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
     1070            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
     1071            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
     1072         !                              ! ---------------------- ! 
     1073      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
     1074         !                              ! ---------------------- ! 
     1075         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
     1076            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
     1077            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
     1078            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
     1079         !                              ! ----------------------- ! 
     1080      ELSE                              ! input cat /= output cat ! 
     1081         !                              ! ----------------------- ! 
     1082          
     1083         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays 
     1084         ALLOCATE( jlmin(idim), jlmax(idim) ) 
     1085 
     1086         ! --- initialize output fields to 0 --- ! 
     1087         ph_i(1:idim,1:jpl) = 0._wp 
     1088         ph_s(1:idim,1:jpl) = 0._wp 
     1089         pa_i(1:idim,1:jpl) = 0._wp 
     1090         ! 
     1091         ! --- fill the categories --- ! 
     1092         !     find where cat-input = cat-output and fill cat-output fields   
     1093         jlmax(:) = 0 
     1094         jlmin(:) = 999 
     1095         jlfil(:,:) = 0 
     1096         DO jl1 = 1, jpl 
     1097            DO jl2 = 1, icat 
     1098               DO ji = 1, idim 
     1099                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 
     1100                     ! fill the right category 
     1101                     ph_i(ji,jl1) = phti(ji,jl2) 
     1102                     ph_s(ji,jl1) = phts(ji,jl2) 
     1103                     pa_i(ji,jl1) = pati(ji,jl2) 
     1104                     ! record categories that are filled 
     1105                     jlmax(ji) = MAX( jlmax(ji), jl1 ) 
     1106                     jlmin(ji) = MIN( jlmin(ji), jl1 ) 
     1107                     jlfil(ji,jl1) = jl1 
     1108                  ENDIF 
     1109               END DO 
     1110            END DO 
     1111         END DO 
     1112         ! 
     1113         ! --- fill the gaps between categories --- !   
     1114         !     transfer from categories filled at the previous step to the empty ones in between 
     1115         DO ji = 1, idim 
     1116            jl1 = jlmin(ji) 
     1117            jl2 = jlmax(ji) 
     1118            IF( jl1 > 1 ) THEN 
     1119               ! fill the lower cat (jl1-1) 
     1120               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 
     1121               ph_i(ji,jl1-1) = hi_mean(jl1-1) 
     1122               ! remove from cat jl1 
     1123               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 
     1124            ENDIF 
     1125            IF( jl2 < jpl ) THEN 
     1126               ! fill the upper cat (jl2+1) 
     1127               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 
     1128               ph_i(ji,jl2+1) = hi_mean(jl2+1) 
     1129               ! remove from cat jl2 
     1130               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 
     1131            ENDIF 
     1132         END DO 
     1133         ! 
     1134         jlfil2(:,:) = jlfil(:,:)  
     1135         ! fill categories from low to high 
     1136         DO jl = 2, jpl-1 
     1137            DO ji = 1, idim 
     1138               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
     1139                  ! fill high 
     1140                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 
     1141                  ph_i(ji,jl) = hi_mean(jl) 
     1142                  jlfil(ji,jl) = jl 
     1143                  ! remove low 
     1144                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 
     1145               ENDIF 
     1146            END DO 
     1147         END DO 
     1148         ! 
     1149         ! fill categories from high to low 
     1150         DO jl = jpl-1, 2, -1 
     1151            DO ji = 1, idim 
     1152               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
     1153                  ! fill low 
     1154                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 
     1155                  ph_i(ji,jl) = hi_mean(jl)  
     1156                  jlfil2(ji,jl) = jl 
     1157                  ! remove high 
     1158                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 
     1159               ENDIF 
     1160            END DO 
     1161         END DO 
     1162         ! 
     1163         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
     1164         DEALLOCATE( jlmin, jlmax ) 
     1165         ! 
     1166         ! == temperature and salinity == ! 
     1167         ! 
     1168         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
     1169         ! 
     1170         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
     1171         ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1172         END WHERE 
     1173         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
     1174         ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
     1175         END WHERE 
     1176         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
     1177         ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
     1178         END WHERE 
     1179         ! 
     1180         ! fill all the categories with the same value 
     1181         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1182         DO jl = 1, jpl 
     1183            pt_i (:,jl) = ztmp(:) 
     1184         END DO 
     1185         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     1186         DO jl = 1, jpl 
     1187            pt_s (:,jl) = ztmp(:) 
     1188         END DO 
     1189         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     1190         DO jl = 1, jpl 
     1191            pt_su(:,jl) = ztmp(:) 
     1192         END DO 
     1193         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1194         DO jl = 1, jpl 
     1195            ps_i (:,jl) = ztmp(:) 
     1196         END DO 
     1197         ! 
     1198         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
     1199         ! 
     1200         ! == ponds == ! 
     1201         ALLOCATE( zfra(idim) ) 
     1202         ! keep the same pond fraction atip/ati for each category 
     1203         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 
     1204         ELSEWHERE                                   ;   zfra(:) = 0._wp 
     1205         END WHERE 
     1206         DO jl = 1, jpl 
     1207            pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1208         END DO 
     1209         ! keep the same v_ip/v_i ratio for each category 
     1210         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1211            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1212         ELSEWHERE 
     1213            zfra(:) = 0._wp 
     1214         END WHERE 
     1215         DO jl = 1, jpl 
     1216            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1217            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     1218            END WHERE 
     1219         END DO 
     1220         DEALLOCATE( zfra ) 
     1221         ! 
     1222      ENDIF 
     1223      ! 
     1224   END SUBROUTINE ice_var_itd_NcMc 
    10401225 
    10411226#else 
Note: See TracChangeset for help on using the changeset viewer.