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 13899 for NEMO/branches/2020/tickets_icb_1900/src/ICE/icevar.F90 – NEMO

Ignore:
Timestamp:
2020-11-27T17:26:33+01:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: update branch to trunk and add ICB test case

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/tickets_icb_1900

    • Property svn:externals
      •  

        old new  
        22^/utils/build/makenemo@HEAD   makenemo 
        33^/utils/build/mk@HEAD         mk 
        4 ^/utils/tools/@HEAD           tools 
         4^/utils/tools@HEAD            tools 
        55^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
         
        88 
        99# SETTE 
        10 ^/utils/CI/sette@12931        sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/tickets_icb_1900/src/ICE/icevar.F90

    r13226 r13899  
    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(:,:,:) 
     
    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         &                                                       ( h_il(:,:,:) - zhl_min ) / ( 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)      
     
    243271      zlay_i   = REAL( nlay_i , wp )    ! number of layers 
    244272      DO jl = 1, jpl 
    245          DO_3D_11_11( 1, nlay_i ) 
     273         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    246274            IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    247275               ! 
     
    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 
     
    347376         z1_dS = 1._wp / ( zsi1 - zsi0 ) 
    348377         DO jl = 1, jpl 
    349             DO_2D_11_11 
     378            DO_2D( 1, 1, 1, 1 ) 
    350379               zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  ) 
    351380               !                             ! force a constant profile when SSS too low (Baltic Sea) 
     
    356385         ! Computation of the profile 
    357386         DO jl = 1, jpl 
    358             DO_3D_11_11( 1, nlay_i ) 
     387            DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    359388               !                          ! linear profile with 0 surface value 
    360389               zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 
     
    486515         ! Zap ice energy and use ocean heat to melt ice 
    487516         !----------------------------------------------------------------- 
    488          DO_3D_11_11( 1, nlay_i ) 
     517         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    489518            ! update exchanges with ocean 
    490519            hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 
     
    493522         END_3D 
    494523         ! 
    495          DO_3D_11_11( 1, nlay_s ) 
     524         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
    496525            ! update exchanges with ocean 
    497526            hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 
     
    503532         ! zap ice and snow volume, add water and salt to ocean 
    504533         !----------------------------------------------------------------- 
    505          DO_2D_11_11 
     534         DO_2D( 1, 1, 1, 1 ) 
    506535            ! update exchanges with ocean 
    507536            sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
     
    521550            a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    522551            v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     552            v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 
    523553            ! 
    524554         END_2D 
     
    542572 
    543573 
    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 ) 
     574   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 ) 
    545575      !!------------------------------------------------------------------- 
    546576      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    557587      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    558588      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     589      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    559590      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    560591      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    574605         ! zap ice energy and send it to the ocean 
    575606         !---------------------------------------- 
    576          DO_3D_11_11( 1, nlay_i ) 
     607         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    577608            IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    578609               hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 
     
    581612         END_3D 
    582613         ! 
    583          DO_3D_11_11( 1, nlay_s ) 
     614         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
    584615            IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    585616               hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 
     
    591622         ! zap ice and snow volume, add water and salt to ocean 
    592623         !----------------------------------------------------- 
    593          DO_2D_11_11 
     624         DO_2D( 1, 1, 1, 1 ) 
    594625            IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
    595626               wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 
     
    613644      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp 
    614645      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 
     646      WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok 
    616647      ! 
    617648   END SUBROUTINE ice_var_zapneg 
    618649 
    619650 
    620    SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     651   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    621652      !!------------------------------------------------------------------- 
    622653      !!                   ***  ROUTINE ice_var_roundoff *** 
     
    631662      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    632663      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     664      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    633665      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    634666      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    643675      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    644676      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    645       IF( ln_pnd_H12 ) THEN 
     677      IF( ln_pnd_LEV ) THEN 
    646678         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    647679         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     680         IF( ln_pnd_lids ) THEN 
     681            WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0 
     682         ENDIF 
    648683      ENDIF 
    649684      ! 
     
    764799   !! ** Purpose :  converting N-cat ice to jpl ice categories 
    765800   !!------------------------------------------------------------------- 
    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 ) 
     801   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     802      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    768803      !!------------------------------------------------------------------- 
    769804      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     
    771806      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    772807      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 
     808      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     809      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 
    775810      !!------------------------------------------------------------------- 
    776811      ! == thickness and concentration == ! 
     
    786821      pa_ip(:) = patip(:) 
    787822      ph_ip(:) = phtip(:) 
     823      ph_il(:) = phtil(:) 
    788824       
    789825   END SUBROUTINE ice_var_itd_1c1c 
    790826 
    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 ) 
     827   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     828      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    793829      !!------------------------------------------------------------------- 
    794830      !! ** Purpose :  converting N-cat ice to 1 ice category 
     
    796832      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    797833      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 
     834      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     835      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 
    800836      ! 
    801837      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
     
    832868      ! == ponds == ! 
    833869      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 
     870      WHERE( pa_ip(:) /= 0._wp ) 
     871         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     872         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     873      ELSEWHERE 
     874         ph_ip(:) = 0._wp 
     875         ph_il(:) = 0._wp 
    836876      END WHERE 
    837877      ! 
     
    840880   END SUBROUTINE ice_var_itd_Nc1c 
    841881    
    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 ) 
     882   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     883      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    844884      !!------------------------------------------------------------------- 
    845885      !! 
     
    863903      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    864904      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 
     905      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     906      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 
    867907      ! 
    868908      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti 
     
    954994         pt_su(:,jl) = ptmsu(:) 
    955995         ps_i (:,jl) = psmi (:) 
    956          ps_i (:,jl) = psmi (:)          
    957996      END DO 
    958997      ! 
     
    9751014         END WHERE 
    9761015      END DO 
     1016      ! keep the same v_il/v_i ratio for each category 
     1017      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     1018      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     1019      END WHERE 
     1020      DO jl = 1, jpl 
     1021         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1022         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp 
     1023         END WHERE 
     1024      END DO 
    9771025      DEALLOCATE( zfra ) 
    9781026      ! 
    9791027   END SUBROUTINE ice_var_itd_1cMc 
    9801028 
    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 ) 
     1029   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, & 
     1030      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 
    9831031      !!------------------------------------------------------------------- 
    9841032      !! 
     
    9951043      !! 
    9961044      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    997        !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     1045      !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    9981046      !!               
    9991047      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
     
    10111059      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    10121060      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 
     1061      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds 
     1062      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 
    10151063      ! 
    10161064      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
     
    10411089         pa_ip(:,:) = patip(:,:) 
    10421090         ph_ip(:,:) = phtip(:,:) 
     1091         ph_il(:,:) = phtil(:,:) 
    10431092         !                              ! ---------------------- ! 
    10441093      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
     
    10461095         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
    10471096            &                    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(:,:)  ) 
     1097            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), & 
     1098            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  ) 
    10501099         !                              ! ---------------------- ! 
    10511100      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
     
    10531102         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
    10541103            &                    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)  ) 
     1104            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), & 
     1105            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  ) 
    10571106         !                              ! ----------------------- ! 
    10581107      ELSE                              ! input cat /= output cat ! 
     
    11961245            END WHERE 
    11971246         END DO 
     1247         ! keep the same v_il/v_i ratio for each category 
     1248         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1249            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1250         ELSEWHERE 
     1251            zfra(:) = 0._wp 
     1252         END WHERE 
     1253         DO jl = 1, jpl 
     1254            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1255            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp 
     1256            END WHERE 
     1257         END DO 
    11981258         DEALLOCATE( zfra ) 
    11991259         ! 
     
    12011261      ! 
    12021262   END SUBROUTINE ice_var_itd_NcMc 
     1263 
     1264   !!------------------------------------------------------------------- 
     1265   !! INTERFACE ice_var_snwfra 
     1266   !! 
     1267   !! ** Purpose :  fraction of ice covered by snow 
     1268   !! 
     1269   !! ** Method  :  In absence of proper snow model on top of sea ice, 
     1270   !!               we argue that snow does not cover the whole ice because 
     1271   !!               of wind blowing... 
     1272   !!                 
     1273   !! ** Arguments : ph_s: snow thickness 
     1274   !!                 
     1275   !! ** Output    : pa_s_fra: fraction of ice covered by snow 
     1276   !! 
     1277   !!------------------------------------------------------------------- 
     1278   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra ) 
     1279      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1280      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1281      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1282         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1283         ELSEWHERE             ; pa_s_fra = 0._wp 
     1284         END WHERE 
     1285      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1286         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1287      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1288         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1289      ENDIF 
     1290   END SUBROUTINE ice_var_snwfra_3d 
     1291 
     1292   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra ) 
     1293      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1294      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1295      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1296         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1297         ELSEWHERE             ; pa_s_fra = 0._wp 
     1298         END WHERE 
     1299      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1300         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1301      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1302         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1303      ENDIF 
     1304   END SUBROUTINE ice_var_snwfra_2d 
     1305 
     1306   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra ) 
     1307      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness 
     1308      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
     1309      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1310         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1311         ELSEWHERE             ; pa_s_fra = 0._wp 
     1312         END WHERE 
     1313      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1314         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1315      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1316         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1317      ENDIF 
     1318   END SUBROUTINE ice_var_snwfra_1d 
     1319    
     1320   !!-------------------------------------------------------------------------- 
     1321   !! INTERFACE ice_var_snwblow 
     1322   !! 
     1323   !! ** Purpose :   Compute distribution of precip over the ice 
     1324   !! 
     1325   !!                Snow accumulation in one thermodynamic time step 
     1326   !!                snowfall is partitionned between leads and ice. 
     1327   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads 
     1328   !!                but because of the winds, more snow falls on leads than on sea ice 
     1329   !!                and a greater fraction (1-at_i)^beta of the total mass of snow  
     1330   !!                (beta < 1) falls in leads. 
     1331   !!                In reality, beta depends on wind speed,  
     1332   !!                and should decrease with increasing wind speed but here, it is  
     1333   !!                considered as a constant. an average value is 0.66 
     1334   !!-------------------------------------------------------------------------- 
     1335!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 
     1336   SUBROUTINE ice_var_snwblow_2d( pin, pout ) 
     1337      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b ) 
     1338      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
     1339      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1340   END SUBROUTINE ice_var_snwblow_2d 
     1341 
     1342   SUBROUTINE ice_var_snwblow_1d( pin, pout ) 
     1343      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
     1344      REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
     1345      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1346   END SUBROUTINE ice_var_snwblow_1d 
    12031347 
    12041348#else 
Note: See TracChangeset for help on using the changeset viewer.