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 11536 for NEMO/trunk/src/ICE/icevar.F90 – NEMO

Ignore:
Timestamp:
2019-09-11T15:54:18+02:00 (5 years ago)
Author:
smasson
Message:

trunk: merge dev_r10984_HPC-13 into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icevar.F90

    r11229 r11536  
    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 
     
    104104      ! 
    105105      !                                      ! integrated values 
    106       vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 ) 
    107       vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 ) 
    108       at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 ) 
    109       et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
    110       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 ) 
    111112      ! 
    112113      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
     
    138139         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    139140         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
    140          sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:) 
     141         sm_i (:,:) =      st_i(:,:)                          * z1_vt_i(:,:) 
    141142         ! 
    142143         tm_i(:,:) = 0._wp 
     
    158159            tm_s (:,:) = rt0 
    159160         END WHERE 
    160  
     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         ! 
    161167         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 
     168         ! 
    162169      ENDIF 
    163170      ! 
     
    263270      ! 
    264271      ! integrated values  
    265       vt_i (:,:) = SUM( v_i, dim=3 ) 
    266       vt_s (:,:) = SUM( v_s, dim=3 ) 
    267       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 ) 
    268275      ! 
    269276   END SUBROUTINE ice_var_glo2eqv 
     
    533540 
    534541      ! to be sure that at_i is the sum of a_i(jl) 
    535       at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 
    536       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 
    537550 
    538551      ! open water = 1 if at_i=0 
     
    652665      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    653666      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    654       IF ( ln_pnd_H12 ) THEN 
     667      IF( ln_pnd_H12 ) THEN 
    655668         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    656669         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     
    773786   !! ** Purpose :  converting N-cat ice to jpl ice categories 
    774787   !!------------------------------------------------------------------- 
    775    SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     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 ) 
    776790      !!------------------------------------------------------------------- 
    777791      !! ** Purpose :  converting 1-cat ice to 1 ice category 
    778792      !!------------------------------------------------------------------- 
    779       REAL(wp), DIMENSION(:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    780       REAL(wp), DIMENSION(:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    781       !!------------------------------------------------------------------- 
    782       zh_i(:) = zhti(:) 
    783       zh_s(:) = zhts(:) 
    784       za_i(:) = zati(:) 
     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       
    785811   END SUBROUTINE ice_var_itd_1c1c 
    786812 
    787    SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     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 ) 
    788815      !!------------------------------------------------------------------- 
    789816      !! ** Purpose :  converting N-cat ice to 1 ice category 
    790817      !!------------------------------------------------------------------- 
    791       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    792       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    793       !!------------------------------------------------------------------- 
    794       ! 
    795       za_i(:) = SUM( zati(:,:), dim=2 ) 
    796       ! 
    797       WHERE( za_i(:) /= 0._wp ) 
    798          zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 
    799          zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 
    800       ELSEWHERE 
    801          zh_i(:) = 0._wp 
    802          zh_s(:) = 0._wp 
     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 
    803837      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 ) 
    804861      ! 
    805862   END SUBROUTINE ice_var_itd_Nc1c 
    806863    
    807    SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     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 ) 
    808866      !!------------------------------------------------------------------- 
    809867      !! 
     
    826884      !!               4) Iterate until ok (SUM(itest(:) = 4) 
    827885      !! 
    828       !! ** Arguments : zhti: 1-cat ice thickness 
    829       !!                zhts: 1-cat snow depth 
    830       !!                zati: 1-cat ice concentration 
     886      !! ** Arguments : phti: 1-cat ice thickness 
     887      !!                phts: 1-cat snow depth 
     888      !!                pati: 1-cat ice concentration 
    831889      !! 
    832890      !! ** Output    : jpl-cat  
     
    834892      !!  (Example of application: BDY forcings when input are cell averaged)   
    835893      !!------------------------------------------------------------------- 
    836       INTEGER  :: ji, jk, jl             ! dummy loop indices 
    837       INTEGER  :: idim, i_fill, jl0   
    838       REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    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(4)                  ::   itest 
    842       !!------------------------------------------------------------------- 
    843       ! 
    844       ! ---------------------------------------- 
     894      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     895      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     896      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     897      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     898      ! 
     899      INTEGER , DIMENSION(4) ::   itest 
     900      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra 
     901      INTEGER  ::   ji, jk, jl 
     902      INTEGER  ::   idim, i_fill, jl0   
     903      REAL(wp) ::   zarg, zV, zconv, zdh, zdv 
     904      !!------------------------------------------------------------------- 
     905      ! 
     906      ! == thickness and concentration == ! 
    845907      ! distribution over the jpl ice categories 
    846       ! ---------------------------------------- 
    847       ! a gaussian distribution for ice concentration is used 
    848       ! then we check whether the distribution fullfills 
    849       ! volume and area conservation, positivity and ice categories bounds 
    850       idim = SIZE( zhti , 1 ) 
    851       zh_i(1:idim,1:jpl) = 0._wp 
    852       zh_s(1:idim,1:jpl) = 0._wp 
    853       za_i(1:idim,1:jpl) = 0._wp 
     908      !    a gaussian distribution for ice concentration is used 
     909      !    then we check whether the distribution fullfills 
     910      !    volume and area conservation, positivity and ice categories bounds 
     911      idim = SIZE( phti , 1 ) 
     912      ! 
     913      ph_i(1:idim,1:jpl) = 0._wp 
     914      ph_s(1:idim,1:jpl) = 0._wp 
     915      pa_i(1:idim,1:jpl) = 0._wp 
    854916      ! 
    855917      DO ji = 1, idim 
    856918         ! 
    857          IF( zhti(ji) > 0._wp ) THEN 
     919         IF( phti(ji) > 0._wp ) THEN 
    858920            ! 
    859921            ! find which category (jl0) the input ice thickness falls into 
    860922            jl0 = jpl 
    861923            DO jl = 1, jpl 
    862                IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     924               IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN 
    863925                  jl0 = jl 
    864926                  CYCLE 
     
    872934               i_fill = i_fill - 1 
    873935               ! 
    874                zh_i(ji,1:jpl) = 0._wp 
    875                za_i(ji,1:jpl) = 0._wp 
     936               ph_i(ji,1:jpl) = 0._wp 
     937               pa_i(ji,1:jpl) = 0._wp 
    876938               itest(:)       = 0       
    877939               ! 
    878940               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    879                   zh_i(ji,1) = zhti(ji) 
    880                   za_i (ji,1) = zati (ji) 
     941                  ph_i(ji,1) = phti(ji) 
     942                  pa_i(ji,1) = pati (ji) 
    881943               ELSE                         !-- case ice is thicker: fill categories >1 
    882944                  ! thickness 
    883945                  DO jl = 1, i_fill - 1 
    884                      zh_i(ji,jl) = hi_mean(jl) 
     946                     ph_i(ji,jl) = hi_mean(jl) 
    885947                  END DO 
    886948                  ! 
    887949                  ! concentration 
    888                   za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 
     950                  pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl)) 
    889951                  DO jl = 1, i_fill - 1 
    890952                     IF ( jl /= jl0 ) THEN 
    891                         zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    892                         za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     953                        zarg        = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp ) 
     954                        pa_i(ji,jl) =   pa_i (ji,jl0) * EXP(-zarg**2) 
    893955                     ENDIF 
    894956                  END DO 
    895957                  ! 
    896958                  ! last category 
    897                   za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) ) 
    898                   zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) ) 
    899                   zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )  
     959                  pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) ) 
     960                  zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) ) 
     961                  ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 )  
    900962                  ! 
    901963                  ! correction if concentration of upper cat is greater than lower cat 
     
    903965                  IF ( jl0 /= jpl ) THEN 
    904966                     DO jl = jpl, jl0+1, -1 
    905                         IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 
    906                            zdv = zh_i(ji,jl) * za_i(ji,jl) 
    907                            zh_i(ji,jl    ) = 0._wp 
    908                            za_i (ji,jl    ) = 0._wp 
    909                            za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 
     967                        IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN 
     968                           zdv = ph_i(ji,jl) * pa_i(ji,jl) 
     969                           ph_i(ji,jl    ) = 0._wp 
     970                           pa_i (ji,jl    ) = 0._wp 
     971                           pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 ) 
    910972                        END IF 
    911973                     END DO 
     
    915977               ! 
    916978               ! Compatibility tests 
    917                zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )  
     979               zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) )  
    918980               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation 
    919981               ! 
    920                zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 
     982               zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) ) 
    921983               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation 
    922984               ! 
    923                IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ? 
     985               IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                   ! Test 3: thickness of the last category is in-bounds ? 
    924986               ! 
    925987               itest(4) = 1 
    926988               DO jl = 1, i_fill 
    927                   IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations 
     989                  IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0                                  ! Test 4: positivity of ice concentrations 
    928990               END DO 
    929991               !                                         !---------------------------- 
     
    933995      END DO 
    934996 
    935       ! Add Snow in each category where za_i is not 0 
     997      ! Add Snow in each category where pa_i is not 0 
    936998      DO jl = 1, jpl 
    937999         DO ji = 1, idim 
    938             IF( za_i(ji,jl) > 0._wp ) THEN 
    939                zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     1000            IF( pa_i(ji,jl) > 0._wp ) THEN 
     1001               ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) ) 
    9401002               ! In case snow load is in excess that would lead to transformation from snow to ice 
    9411003               ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    942                zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )  
     1004               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )  
    9431005               ! recompute h_i, h_s avoiding out of bounds values 
    944                zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh ) 
    945                zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos ) 
     1006               ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 
     1007               ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 
    9461008            ENDIF 
    9471009         END DO 
    9481010      END DO 
    9491011      ! 
     1012      ! == temperature and salinity == ! 
     1013      DO jl = 1, jpl 
     1014         pt_i (:,jl) = ptmi (:) 
     1015         pt_s (:,jl) = ptms (:) 
     1016         pt_su(:,jl) = ptmsu(:) 
     1017         ps_i (:,jl) = psmi (:) 
     1018         ps_i (:,jl) = psmi (:)          
     1019      END DO 
     1020      ! 
     1021      ! == ponds == ! 
     1022      ALLOCATE( zfra(idim) ) 
     1023      ! keep the same pond fraction atip/ati for each category 
     1024      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:) 
     1025      ELSEWHERE                   ;   zfra(:) = 0._wp 
     1026      END WHERE 
     1027      DO jl = 1, jpl 
     1028         pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1029      END DO 
     1030      ! keep the same v_ip/v_i ratio for each category 
     1031      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     1032      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     1033      END WHERE 
     1034      DO jl = 1, jpl 
     1035         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1036         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     1037         END WHERE 
     1038      END DO 
     1039      DEALLOCATE( zfra ) 
     1040      ! 
    9501041   END SUBROUTINE ice_var_itd_1cMc 
    9511042 
    952    SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     1043   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     1044      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    9531045      !!------------------------------------------------------------------- 
    9541046      !! 
     
    9711063      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 
    9721064      !! 
    973       !! ** Arguments : zhti: N-cat ice thickness 
    974       !!                zhts: N-cat snow depth 
    975       !!                zati: N-cat ice concentration 
     1065      !! ** Arguments : phti: N-cat ice thickness 
     1066      !!                phts: N-cat snow depth 
     1067      !!                pati: N-cat ice concentration 
    9761068      !! 
    9771069      !! ** Output    : jpl-cat  
     
    9791071      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)   
    9801072      !!------------------------------------------------------------------- 
    981       INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices 
     1073      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     1074      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
     1075      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     1076      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     1077      ! 
     1078      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
     1079      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin 
     1080      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra 
     1081      ! 
     1082      REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
     1083      INTEGER  ::   ji, jl, jl1, jl2 
    9821084      INTEGER  ::   idim, icat   
    983       REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
    984       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    985       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    986       INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2 
    987       INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin 
    988       !!------------------------------------------------------------------- 
    989       ! 
    990       idim = SIZE( zhti, 1 ) 
    991       icat = SIZE( zhti, 2 ) 
     1085      !!------------------------------------------------------------------- 
     1086      ! 
     1087      idim = SIZE( phti, 1 ) 
     1088      icat = SIZE( phti, 2 ) 
     1089      ! 
     1090      ! == thickness and concentration == ! 
    9921091      !                                 ! ---------------------- ! 
    9931092      IF( icat == jpl ) THEN            ! input cat = output cat ! 
    9941093         !                              ! ---------------------- ! 
    995          zh_i(:,:) = zhti(:,:) 
    996          zh_s(:,:) = zhts(:,:) 
    997          za_i(:,:) = zati(:,:) 
     1094         ph_i(:,:) = phti(:,:) 
     1095         ph_s(:,:) = phts(:,:) 
     1096         pa_i(:,:) = pati(:,:) 
     1097         ! 
     1098         ! == temperature and salinity and ponds == ! 
     1099         pt_i (:,:) = ptmi (:,:) 
     1100         pt_s (:,:) = ptms (:,:) 
     1101         pt_su(:,:) = ptmsu(:,:) 
     1102         ps_i (:,:) = psmi (:,:) 
     1103         pa_ip(:,:) = patip(:,:) 
     1104         ph_ip(:,:) = phtip(:,:) 
    9981105         !                              ! ---------------------- ! 
    9991106      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
    10001107         !                              ! ---------------------- ! 
    1001          CALL  ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) ) 
     1108         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
     1109            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
     1110            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
     1111            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
    10021112         !                              ! ---------------------- ! 
    10031113      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
    10041114         !                              ! ---------------------- ! 
    1005          CALL  ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) )          
     1115         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
     1116            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
     1117            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
     1118            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
    10061119         !                              ! ----------------------- ! 
    10071120      ELSE                              ! input cat /= output cat ! 
     
    10121125 
    10131126         ! --- initialize output fields to 0 --- ! 
    1014          zh_i(1:idim,1:jpl) = 0._wp 
    1015          zh_s(1:idim,1:jpl) = 0._wp 
    1016          za_i(1:idim,1:jpl) = 0._wp 
     1127         ph_i(1:idim,1:jpl) = 0._wp 
     1128         ph_s(1:idim,1:jpl) = 0._wp 
     1129         pa_i(1:idim,1:jpl) = 0._wp 
    10171130         ! 
    10181131         ! --- fill the categories --- ! 
     
    10241137            DO jl2 = 1, icat 
    10251138               DO ji = 1, idim 
    1026                   IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
     1139                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 
    10271140                     ! fill the right category 
    1028                      zh_i(ji,jl1) = zhti(ji,jl2) 
    1029                      zh_s(ji,jl1) = zhts(ji,jl2) 
    1030                      za_i(ji,jl1) = zati(ji,jl2) 
     1141                     ph_i(ji,jl1) = phti(ji,jl2) 
     1142                     ph_s(ji,jl1) = phts(ji,jl2) 
     1143                     pa_i(ji,jl1) = pati(ji,jl2) 
    10311144                     ! record categories that are filled 
    10321145                     jlmax(ji) = MAX( jlmax(ji), jl1 ) 
     
    10451158            IF( jl1 > 1 ) THEN 
    10461159               ! fill the lower cat (jl1-1) 
    1047                za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 
    1048                zh_i(ji,jl1-1) = hi_mean(jl1-1) 
     1160               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 
     1161               ph_i(ji,jl1-1) = hi_mean(jl1-1) 
    10491162               ! remove from cat jl1 
    1050                za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
     1163               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 
    10511164            ENDIF 
    10521165            IF( jl2 < jpl ) THEN 
    10531166               ! fill the upper cat (jl2+1) 
    1054                za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 
    1055                zh_i(ji,jl2+1) = hi_mean(jl2+1) 
     1167               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 
     1168               ph_i(ji,jl2+1) = hi_mean(jl2+1) 
    10561169               ! remove from cat jl2 
    1057                za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
     1170               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 
    10581171            ENDIF 
    10591172         END DO 
     
    10651178               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
    10661179                  ! fill high 
    1067                   za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
    1068                   zh_i(ji,jl) = hi_mean(jl) 
     1180                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 
     1181                  ph_i(ji,jl) = hi_mean(jl) 
    10691182                  jlfil(ji,jl) = jl 
    10701183                  ! remove low 
    1071                   za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
     1184                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 
    10721185               ENDIF 
    10731186            END DO 
     
    10791192               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
    10801193                  ! fill low 
    1081                   za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 
    1082                   zh_i(ji,jl) = hi_mean(jl)  
     1194                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 
     1195                  ph_i(ji,jl) = hi_mean(jl)  
    10831196                  jlfil2(ji,jl) = jl 
    10841197                  ! remove high 
    1085                   za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
     1198                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 
    10861199               ENDIF 
    10871200            END DO 
     
    10901203         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
    10911204         DEALLOCATE( jlmin, jlmax ) 
     1205         ! 
     1206         ! == temperature and salinity == ! 
     1207         ! 
     1208         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
     1209         ! 
     1210         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
     1211         ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1212         END WHERE 
     1213         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
     1214         ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
     1215         END WHERE 
     1216         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
     1217         ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
     1218         END WHERE 
     1219         ! 
     1220         ! fill all the categories with the same value 
     1221         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1222         DO jl = 1, jpl 
     1223            pt_i (:,jl) = ztmp(:) 
     1224         END DO 
     1225         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     1226         DO jl = 1, jpl 
     1227            pt_s (:,jl) = ztmp(:) 
     1228         END DO 
     1229         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     1230         DO jl = 1, jpl 
     1231            pt_su(:,jl) = ztmp(:) 
     1232         END DO 
     1233         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1234         DO jl = 1, jpl 
     1235            ps_i (:,jl) = ztmp(:) 
     1236         END DO 
     1237         ! 
     1238         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
     1239         ! 
     1240         ! == ponds == ! 
     1241         ALLOCATE( zfra(idim) ) 
     1242         ! keep the same pond fraction atip/ati for each category 
     1243         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 
     1244         ELSEWHERE                                   ;   zfra(:) = 0._wp 
     1245         END WHERE 
     1246         DO jl = 1, jpl 
     1247            pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1248         END DO 
     1249         ! keep the same v_ip/v_i ratio for each category 
     1250         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1251            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1252         ELSEWHERE 
     1253            zfra(:) = 0._wp 
     1254         END WHERE 
     1255         DO jl = 1, jpl 
     1256            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1257            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     1258            END WHERE 
     1259         END DO 
     1260         DEALLOCATE( zfra ) 
    10921261         ! 
    10931262      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.