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 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icevar.F90 – NEMO

Ignore:
Timestamp:
2019-09-19T11:18:03+02:00 (5 years ago)
Author:
jchanut
Message:

#2222, merged with trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icevar.F90

    r11229 r11573  
    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      !! 
    810868      !! ** Purpose :  converting 1-cat ice to jpl ice categories 
    811869      !! 
    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) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
    819       !! 
    820       !!               2) Check whether the distribution conserves area and volume, positivity and 
    821       !!                  category boundaries 
     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 
    822873      !!               
    823       !!               3) If not (input ice is too thin), the last category is empty and 
    824       !!                  the number of categories is reduced (jpl-1) 
    825       !! 
    826       !!               4) Iterate until ok (SUM(itest(:) = 4) 
    827       !! 
    828       !! ** Arguments : zhti: 1-cat ice thickness 
    829       !!                zhts: 1-cat snow depth 
    830       !!                zati: 1-cat ice concentration 
     874      !! 
     875      !! ** Arguments : phti: 1-cat ice thickness 
     876      !!                phts: 1-cat snow depth 
     877      !!                pati: 1-cat ice concentration 
    831878      !! 
    832879      !! ** Output    : jpl-cat  
    833880      !! 
    834       !!  (Example of application: BDY forcings when input are cell averaged)   
    835       !!------------------------------------------------------------------- 
    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       ! ---------------------------------------- 
    845       ! 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 
    854       ! 
     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 
    855935      DO ji = 1, idim 
    856936         ! 
    857          IF( zhti(ji) > 0._wp ) THEN 
    858             ! 
    859             ! find which category (jl0) the input ice thickness falls into 
    860             jl0 = jpl 
    861             DO jl = 1, jpl 
    862                IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
    863                   jl0 = jl 
    864                   CYCLE 
    865                ENDIF 
    866             END DO 
    867             ! 
    868             itest(:) = 0 
    869             i_fill   = jpl + 1                                            !------------------------------------ 
    870             DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    871                !                                                          !------------------------------------ 
    872                i_fill = i_fill - 1 
    873                ! 
    874                zh_i(ji,1:jpl) = 0._wp 
    875                za_i(ji,1:jpl) = 0._wp 
    876                itest(:)       = 0       
    877                ! 
    878                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) 
    881                ELSE                         !-- case ice is thicker: fill categories >1 
    882                   ! thickness 
    883                   DO jl = 1, i_fill - 1 
    884                      zh_i(ji,jl) = hi_mean(jl) 
    885                   END DO 
    886                   ! 
    887                   ! concentration 
    888                   za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 
    889                   DO jl = 1, i_fill - 1 
    890                      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) 
    893                      ENDIF 
    894                   END DO 
    895                   ! 
    896                   ! 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 )  
    900                   ! 
    901                   ! correction if concentration of upper cat is greater than lower cat 
    902                   !    (it should be a gaussian around jl0 but sometimes it is not) 
    903                   IF ( jl0 /= jpl ) THEN 
    904                      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 ) 
    910                         END IF 
    911                      END DO 
    912                   ENDIF 
    913                   ! 
    914                ENDIF 
    915                ! 
    916                ! Compatibility tests 
    917                zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )  
    918                IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation 
    919                ! 
    920                zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 
    921                IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation 
    922                ! 
    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 ? 
    924                ! 
    925                itest(4) = 1 
    926                DO jl = 1, i_fill 
    927                   IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations 
    928                END DO 
    929                !                                         !---------------------------- 
    930             END DO                                       ! end iteration on categories 
    931             !                                            !---------------------------- 
     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 
    932951         ENDIF 
    933       END DO 
    934  
    935       ! Add Snow in each category where za_i is not 0 
     952         ! 
     953      ENDDO 
     954      ! 
     955      ! Add Snow in each category where pa_i is not 0 
    936956      DO jl = 1, jpl 
    937957         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) ) 
     958            IF( pa_i(ji,jl) > 0._wp ) THEN 
     959               ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 
    940960               ! In case snow load is in excess that would lead to transformation from snow to ice 
    941961               ! 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 )  
     962               zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )  
    943963               ! 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 ) 
     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 ) 
    946966            ENDIF 
    947967         END DO 
    948968      END DO 
    949969      ! 
     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      ! 
    9501001   END SUBROUTINE ice_var_itd_1cMc 
    9511002 
    952    SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     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 ) 
    9531005      !!------------------------------------------------------------------- 
    9541006      !! 
     
    9711023      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 
    9721024      !! 
    973       !! ** Arguments : zhti: N-cat ice thickness 
    974       !!                zhts: N-cat snow depth 
    975       !!                zati: N-cat ice concentration 
     1025      !! ** Arguments : phti: N-cat ice thickness 
     1026      !!                phts: N-cat snow depth 
     1027      !!                pati: N-cat ice concentration 
    9761028      !! 
    9771029      !! ** Output    : jpl-cat  
     
    9791031      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)   
    9801032      !!------------------------------------------------------------------- 
    981       INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices 
     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 
    9821044      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 ) 
     1045      !!------------------------------------------------------------------- 
     1046      ! 
     1047      idim = SIZE( phti, 1 ) 
     1048      icat = SIZE( phti, 2 ) 
     1049      ! 
     1050      ! == thickness and concentration == ! 
    9921051      !                                 ! ---------------------- ! 
    9931052      IF( icat == jpl ) THEN            ! input cat = output cat ! 
    9941053         !                              ! ---------------------- ! 
    995          zh_i(:,:) = zhti(:,:) 
    996          zh_s(:,:) = zhts(:,:) 
    997          za_i(:,:) = zati(:,:) 
     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(:,:) 
    9981065         !                              ! ---------------------- ! 
    9991066      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
    10001067         !                              ! ---------------------- ! 
    1001          CALL  ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) ) 
     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(:,:)  ) 
    10021072         !                              ! ---------------------- ! 
    10031073      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
    10041074         !                              ! ---------------------- ! 
    1005          CALL  ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) )          
     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)  ) 
    10061079         !                              ! ----------------------- ! 
    10071080      ELSE                              ! input cat /= output cat ! 
     
    10121085 
    10131086         ! --- 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 
     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 
    10171090         ! 
    10181091         ! --- fill the categories --- ! 
     
    10241097            DO jl2 = 1, icat 
    10251098               DO ji = 1, idim 
    1026                   IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
     1099                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 
    10271100                     ! 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) 
     1101                     ph_i(ji,jl1) = phti(ji,jl2) 
     1102                     ph_s(ji,jl1) = phts(ji,jl2) 
     1103                     pa_i(ji,jl1) = pati(ji,jl2) 
    10311104                     ! record categories that are filled 
    10321105                     jlmax(ji) = MAX( jlmax(ji), jl1 ) 
     
    10451118            IF( jl1 > 1 ) THEN 
    10461119               ! 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) 
     1120               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 
     1121               ph_i(ji,jl1-1) = hi_mean(jl1-1) 
    10491122               ! remove from cat jl1 
    1050                za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
     1123               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 
    10511124            ENDIF 
    10521125            IF( jl2 < jpl ) THEN 
    10531126               ! 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) 
     1127               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 
     1128               ph_i(ji,jl2+1) = hi_mean(jl2+1) 
    10561129               ! remove from cat jl2 
    1057                za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
     1130               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 
    10581131            ENDIF 
    10591132         END DO 
     
    10651138               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
    10661139                  ! fill high 
    1067                   za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
    1068                   zh_i(ji,jl) = hi_mean(jl) 
     1140                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 
     1141                  ph_i(ji,jl) = hi_mean(jl) 
    10691142                  jlfil(ji,jl) = jl 
    10701143                  ! remove low 
    1071                   za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
     1144                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 
    10721145               ENDIF 
    10731146            END DO 
     
    10791152               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
    10801153                  ! 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)  
     1154                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 
     1155                  ph_i(ji,jl) = hi_mean(jl)  
    10831156                  jlfil2(ji,jl) = jl 
    10841157                  ! remove high 
    1085                   za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
     1158                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 
    10861159               ENDIF 
    10871160            END DO 
     
    10901163         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
    10911164         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 ) 
    10921221         ! 
    10931222      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.