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 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icevar.F90 – NEMO

Ignore:
Timestamp:
2020-12-03T12:20:38+01:00 (3 years ago)
Author:
ayoung
Message:

Updated to trunk at 14020. Sette tests passed with change of results for configurations with non-linear ssh. Ticket #2506.

Location:
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13292        sette 
         10^/utils/CI/sette_wave@13990         sette 
  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icevar.F90

    r13295 r14037  
    5151   !!   ice_var_sshdyn    : compute equivalent ssh in lead 
    5252   !!   ice_var_itd       : convert N-cat to M-cat 
     53   !!   ice_var_snwfra    : fraction of ice covered by snow 
     54   !!   ice_var_snwblow   : distribute snow fall between ice and ocean 
    5355   !!---------------------------------------------------------------------- 
    5456   USE dom_oce        ! ocean space and time domain 
     
    7779   PUBLIC   ice_var_sshdyn 
    7880   PUBLIC   ice_var_itd 
     81   PUBLIC   ice_var_snwfra 
     82   PUBLIC   ice_var_snwblow 
    7983 
    8084   INTERFACE ice_var_itd 
     
    8488   !! * Substitutions 
    8589#  include "do_loop_substitute.h90" 
     90 
     91   INTERFACE ice_var_snwfra 
     92      MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d 
     93   END INTERFACE 
     94 
     95   INTERFACE ice_var_snwblow 
     96      MODULE PROCEDURE ice_var_snwblow_1d, ice_var_snwblow_2d 
     97   END INTERFACE 
     98 
    8699   !!---------------------------------------------------------------------- 
    87100   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    115128      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
    116129      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
     130      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) 
    117131      ! 
    118132      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction   
     
    166180         ! 
    167181         !                           ! mean melt pond depth 
    168          WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 
    169          ELSEWHERE                      ;   hm_ip(:,:) = 0._wp 
     182         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)   ;   hm_il(:,:) = vt_il(:,:) / at_ip(:,:) 
     183         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp                     ;   hm_il(:,:) = 0._wp 
    170184         END WHERE          
    171185         ! 
     
    191205      REAL(wp) ::   zhmax, z1_zhmax                 !   -      - 
    192206      REAL(wp) ::   zlay_i, zlay_s                  !   -      - 
    193       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i 
     207      REAL(wp), PARAMETER ::   zhl_max =  0.015_wp  ! pond lid thickness above which the ponds disappear from the albedo calculation 
     208      REAL(wp), PARAMETER ::   zhl_min =  0.005_wp  ! pond lid thickness below which the full pond area is used in the albedo calculation 
     209      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i, z1_a_ip, za_s_fra 
    194210      !!------------------------------------------------------------------- 
    195211 
     
    210226      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp 
    211227      END WHERE 
     228      ! 
     229      WHERE( a_ip(:,:,:) > epsi20 )  ;   z1_a_ip(:,:,:) = 1._wp / a_ip(:,:,:) 
     230      ELSEWHERE                      ;   z1_a_ip(:,:,:) = 0._wp 
     231      END WHERE 
    212232      !                                           !--- ice thickness 
    213233      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 
     
    216236      z1_zhmax =  1._wp / hi_max(jpl)                
    217237      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
    218          h_i  (:,:,jpl) = zhmax 
     238         h_i   (:,:,jpl) = zhmax 
    219239         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    220240         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 
     
    224244      !                                           !--- ice age       
    225245      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) 
    226       !                                           !--- pond fraction and thickness       
     246      !                                           !--- pond and lid thickness       
     247      h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:) 
     248      h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:) 
     249      !                                           !--- melt pond effective area (used for albedo) 
    227250      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:) 
    228       WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:) 
    229       ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp 
    230       END WHERE 
     251      WHERE    ( h_il(:,:,:) <= zhl_min )  ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:)       ! lid is very thin.  Expose all the pond 
     252      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
     253      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond 
     254         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
     255      END WHERE 
     256      ! 
     257      CALL ice_var_snwfra( h_s, za_s_fra )           ! calculate ice fraction covered by snow 
     258      a_ip_eff = MIN( a_ip_eff, 1._wp - za_s_fra )   ! make sure (a_ip_eff + a_s_fra) <= 1 
    231259      ! 
    232260      !                                           !---  salinity (with a minimum value imposed everywhere)      
     
    292320      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:) 
    293321      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
     322      v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 
    294323      ! 
    295324   END SUBROUTINE ice_var_eqv2glo 
     
    505534         DO_2D( 1, 1, 1, 1 ) 
    506535            ! update exchanges with ocean 
    507             sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
    508             wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
    509             wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
     536            sfx_res(ji,jj)  = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
     537            wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
     538            wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
     539            wfx_pnd(ji,jj)  = wfx_pnd(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 
    510540            ! 
    511541            a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 
     
    521551            a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    522552            v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     553            v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 
     554            h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 
     555            h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 
    523556            ! 
    524557         END_2D 
     
    542575 
    543576 
    544    SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     577   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    545578      !!------------------------------------------------------------------- 
    546579      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    557590      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    558591      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     592      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    559593      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    560594      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    604638               psv_i  (ji,jj,jl) = 0._wp 
    605639            ENDIF 
     640            IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 
     641               wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 
     642               pv_il  (ji,jj,jl) = 0._wp 
     643            ENDIF 
     644            IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 
     645               wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 
     646               pv_ip  (ji,jj,jl) = 0._wp 
     647            ENDIF 
    606648         END_2D 
    607649         ! 
     
    612654      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp 
    613655      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp 
    614       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    615       !                                                        but it does not change conservation, so keep it this way is ok 
    616656      ! 
    617657   END SUBROUTINE ice_var_zapneg 
    618658 
    619659 
    620    SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     660   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    621661      !!------------------------------------------------------------------- 
    622662      !!                   ***  ROUTINE ice_var_roundoff *** 
     
    631671      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    632672      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     673      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    633674      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    634675      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    643684      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    644685      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    645       IF( ln_pnd_H12 ) THEN 
     686      IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    646687         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    647688         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     689         IF( ln_pnd_lids ) THEN 
     690            WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0 
     691         ENDIF 
    648692      ENDIF 
    649693      ! 
     
    764808   !! ** Purpose :  converting N-cat ice to jpl ice categories 
    765809   !!------------------------------------------------------------------- 
    766    SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    767       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     810   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     811      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    768812      !!------------------------------------------------------------------- 
    769813      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     
    771815      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    772816      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    773       REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    774       REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     817      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     818      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds 
    775819      !!------------------------------------------------------------------- 
    776820      ! == thickness and concentration == ! 
     
    786830      pa_ip(:) = patip(:) 
    787831      ph_ip(:) = phtip(:) 
     832      ph_il(:) = phtil(:) 
    788833       
    789834   END SUBROUTINE ice_var_itd_1c1c 
    790835 
    791    SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    792       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     836   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     837      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    793838      !!------------------------------------------------------------------- 
    794839      !! ** Purpose :  converting N-cat ice to 1 ice category 
     
    796841      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    797842      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    798       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    799       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     843      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     844      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds 
    800845      ! 
    801846      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
     
    832877      ! == ponds == ! 
    833878      pa_ip(:) = SUM( patip(:,:), dim=2 ) 
    834       WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
    835       ELSEWHERE                    ;   ph_ip(:) = 0._wp 
     879      WHERE( pa_ip(:) /= 0._wp ) 
     880         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     881         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     882      ELSEWHERE 
     883         ph_ip(:) = 0._wp 
     884         ph_il(:) = 0._wp 
    836885      END WHERE 
    837886      ! 
     
    840889   END SUBROUTINE ice_var_itd_Nc1c 
    841890    
    842    SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    843       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     891   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     892      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    844893      !!------------------------------------------------------------------- 
    845894      !! 
     
    863912      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    864913      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    865       REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    866       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     914      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     915      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds 
    867916      ! 
    868917      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti 
     
    9541003         pt_su(:,jl) = ptmsu(:) 
    9551004         ps_i (:,jl) = psmi (:) 
    956          ps_i (:,jl) = psmi (:)          
    9571005      END DO 
    9581006      ! 
     
    9751023         END WHERE 
    9761024      END DO 
     1025      ! keep the same v_il/v_i ratio for each category 
     1026      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     1027      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     1028      END WHERE 
     1029      DO jl = 1, jpl 
     1030         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1031         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp 
     1032         END WHERE 
     1033      END DO 
    9771034      DEALLOCATE( zfra ) 
    9781035      ! 
    9791036   END SUBROUTINE ice_var_itd_1cMc 
    9801037 
    981    SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
    982       &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
     1038   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     1039      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    9831040      !!------------------------------------------------------------------- 
    9841041      !! 
     
    9951052      !! 
    9961053      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    997        !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     1054      !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    9981055      !!               
    9991056      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
     
    10111068      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    10121069      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    1013       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
    1014       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
     1070      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     1071      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds 
    10151072      ! 
    10161073      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
     
    10411098         pa_ip(:,:) = patip(:,:) 
    10421099         ph_ip(:,:) = phtip(:,:) 
     1100         ph_il(:,:) = phtil(:,:) 
    10431101         !                              ! ---------------------- ! 
    10441102      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
     
    10461104         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
    10471105            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
    1048             &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
    1049             &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
     1106            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), & 
     1107            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  ) 
    10501108         !                              ! ---------------------- ! 
    10511109      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
     
    10531111         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
    10541112            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
    1055             &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
    1056             &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
     1113            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), & 
     1114            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  ) 
    10571115         !                              ! ----------------------- ! 
    10581116      ELSE                              ! input cat /= output cat ! 
     
    11961254            END WHERE 
    11971255         END DO 
     1256         ! keep the same v_il/v_i ratio for each category 
     1257         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1258            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1259         ELSEWHERE 
     1260            zfra(:) = 0._wp 
     1261         END WHERE 
     1262         DO jl = 1, jpl 
     1263            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1264            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp 
     1265            END WHERE 
     1266         END DO 
    11981267         DEALLOCATE( zfra ) 
    11991268         ! 
     
    12011270      ! 
    12021271   END SUBROUTINE ice_var_itd_NcMc 
     1272 
     1273   !!------------------------------------------------------------------- 
     1274   !! INTERFACE ice_var_snwfra 
     1275   !! 
     1276   !! ** Purpose :  fraction of ice covered by snow 
     1277   !! 
     1278   !! ** Method  :  In absence of proper snow model on top of sea ice, 
     1279   !!               we argue that snow does not cover the whole ice because 
     1280   !!               of wind blowing... 
     1281   !!                 
     1282   !! ** Arguments : ph_s: snow thickness 
     1283   !!                 
     1284   !! ** Output    : pa_s_fra: fraction of ice covered by snow 
     1285   !! 
     1286   !!------------------------------------------------------------------- 
     1287   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra ) 
     1288      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1289      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1290      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1291         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1292         ELSEWHERE             ; pa_s_fra = 0._wp 
     1293         END WHERE 
     1294      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1295         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1296      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1297         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1298      ENDIF 
     1299   END SUBROUTINE ice_var_snwfra_3d 
     1300 
     1301   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra ) 
     1302      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1303      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1304      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1305         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1306         ELSEWHERE             ; pa_s_fra = 0._wp 
     1307         END WHERE 
     1308      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1309         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1310      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1311         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1312      ENDIF 
     1313   END SUBROUTINE ice_var_snwfra_2d 
     1314 
     1315   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra ) 
     1316      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1317      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1318      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1319         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1320         ELSEWHERE             ; pa_s_fra = 0._wp 
     1321         END WHERE 
     1322      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1323         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1324      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1325         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1326      ENDIF 
     1327   END SUBROUTINE ice_var_snwfra_1d 
     1328    
     1329   !!-------------------------------------------------------------------------- 
     1330   !! INTERFACE ice_var_snwblow 
     1331   !! 
     1332   !! ** Purpose :   Compute distribution of precip over the ice 
     1333   !! 
     1334   !!                Snow accumulation in one thermodynamic time step 
     1335   !!                snowfall is partitionned between leads and ice. 
     1336   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads 
     1337   !!                but because of the winds, more snow falls on leads than on sea ice 
     1338   !!                and a greater fraction (1-at_i)^beta of the total mass of snow  
     1339   !!                (beta < 1) falls in leads. 
     1340   !!                In reality, beta depends on wind speed,  
     1341   !!                and should decrease with increasing wind speed but here, it is  
     1342   !!                considered as a constant. an average value is 0.66 
     1343   !!-------------------------------------------------------------------------- 
     1344!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 
     1345   SUBROUTINE ice_var_snwblow_2d( pin, pout ) 
     1346      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b ) 
     1347      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
     1348      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1349   END SUBROUTINE ice_var_snwblow_2d 
     1350 
     1351   SUBROUTINE ice_var_snwblow_1d( pin, pout ) 
     1352      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
     1353      REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
     1354      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1355   END SUBROUTINE ice_var_snwblow_1d 
    12031356 
    12041357#else 
Note: See TracChangeset for help on using the changeset viewer.