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 14005 – NEMO

Changeset 14005


Ignore:
Timestamp:
2020-12-02T15:16:36+01:00 (3 years ago)
Author:
clem
Message:

merge branch SI3_martin_ponds into the trunk

Location:
NEMO/trunk
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/cfgs/SHARED/field_def_nemo-ice.xml

    r13999 r14005  
    5151          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" />  
    5252          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" />  
     53          <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="kg/m2/s" /> 
     54          <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="kg/m2/s" /> 
     55          <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="kg/m2/s" /> 
     56          <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="kg/m2/s" />   
    5357      
    5458     <!-- heat --> 
    … …  
    302306          <field id="snwtemp_cat"  long_name="Snow temperature per category"                     unit="degC"    detect_missing_value="true" /> 
    303307          <field id="icettop_cat"  long_name="Ice/snow surface temperature per category"         unit="degC"    detect_missing_value="true" /> 
    304           <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        />  
     308          <field id="iceapnd_cat"  long_name="Ice melt pond grid fraction per category"          unit=""        />  
     309          <field id="icevpnd_cat"  long_name="Ice melt pond volume per grid area per category"   unit="m"       />   
    305310          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
    306311          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    307           <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     312          <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category"           unit=""        />  
    308313          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
    309314          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 

  • NEMO/trunk/cfgs/SHARED/namelist_ice_ref

    r13999 r14005  
    2424   jpl              =   5             !  number of ice  categories 
    2525   nlay_i           =   2             !  number of ice  layers 
    26    nlay_s           =   1             !  number of snow layers (only 1 is working) 
     26   nlay_s           =   2             !  number of snow layers 
    2727   ln_virtual_itd   =   .false.       !  virtual ITD mono-category parameterization (jpl=1 only) 
    2828                                      !     i.e. enhanced thermal conductivity & virtual thin ice melting 
    … …  
    6262      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
    6363      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??] 
     64 
     65   cn_dir           = './'      !  root directory for the grounded icebergs mask data location 
     66   !___________!________________!___________________!___________!_____________!________!___________!__________!__________!_______________! 
     67   !           !  file name     ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     68   !           !                !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' ! filename ! pairing  !    filename   ! 
     69   sn_icbmsk       = 'NOT USED' ,       -12.        , 'icb_mask',   .false.   , .true. , 'yearly'  , ''       , ''       , '' 
    6470/ 
    6571!------------------------------------------------------------------------------ 
    … …  
    196202!------------------------------------------------------------------------------ 
    197203   ln_pnd            = .true.         !  activate melt ponds or not 
    198       ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 
    199          rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 
    200          rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 
     204      ln_pnd_TOPO    = .false.        !  topographic melt ponds 
     205      ln_pnd_LEV     = .true.         !  level ice melt ponds 
     206         rn_apnd_min =   0.15         !     minimum meltwater fraction contributing to pond growth (TOPO and LEV) 
     207         rn_apnd_max =   0.85         !     maximum meltwater fraction contributing to pond growth (TOPO and LEV) 
     208         rn_pnd_flush=   0.01         !     pond flushing efficiency (tuning parameter) (LEV) 
    201209      ln_pnd_CST     = .false.        !  constant  melt ponds 
    202210         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC 
    … …  
    262270   ln_icediachk     = .false.         !  check online heat, mass & salt budgets 
    263271      !                               !   rate of ice spuriously gained/lost at each time step => rn_icechk=1 <=> 1.e-6 m/hour 
    264       rn_icechk_cel =  100.           !     check at each gridcell          (1.e-4m/h)=> stops the code if violated (and writes a file) 
    265       rn_icechk_glo =  1.             !     check over the entire ice cover (1.e-6m/h)=> only prints warnings 
     272      rn_icechk_cel =  1.             !     check at each gridcell          (1.e-06m/h)=> stops the code if violated (and writes a file) 
     273      rn_icechk_glo =  1.e-04         !     check over the entire ice cover (1.e-10m/h)=> only prints warnings 
    266274   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    267275   ln_icectl        = .false.         !  ice points output for debug (T or F) 

  • NEMO/trunk/src/ICE/ice.F90

    r13999 r14005  
    210210   !                                     !!** ice-ponds namelist (namthd_pnd) 
    211211   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    212    LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 
    213    REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds 
    214    REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds 
     212   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 
     213   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Simple melt pond scheme 
     214   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum fraction of melt water contributing to ponds 
     215   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum fraction of melt water contributing to ponds 
     216   REAL(wp), PUBLIC ::   rn_pnd_flush     !: Pond flushing efficiency (tuning parameter) 
    215217   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    216218   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
    … …  
    345347   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   om_i          !: mean ice age over all categories                        (s) 
    346348   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tau_icebfr    !: ice friction on ocean bottom (landfast param activated) 
     349   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   icb_mask      !: mask of grounded icebergs if landfast [0-1] 
    347350 
    348351   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s           !: Snow temperatures     [K] 
    … …  
    366369   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   vt_il         !: total melt pond lid volume per gridcell area [m] 
    367370 
     371   ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array) 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
     373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
     374 
    368375   !!---------------------------------------------------------------------- 
    369376   !! * Global variables at before time step 
    370377   !!---------------------------------------------------------------------- 
    371378   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 
     379   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_ip_b, v_il_b             !: ponds and lids volumes 
    372380   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b              !: 
    373381   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                      !: snow heat content 
    … …  
    396404   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
    397405   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1]  
     406   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vpnd         !: pond volume variation  [m/s]  
    398407   ! 
    399408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s) 
    … …  
    473482         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) ,  & 
    474483         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) ,  & 
    475          &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj)            , STAT=ierr(ii) ) 
     484         &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), icb_mask(jpi,jpj), STAT=ierr(ii) ) 
    476485 
    477486      ii = ii + 1 
    … …  
    483492      ii = ii + 1 
    484493      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  & 
    485          &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 
     494         &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) ,                     & 
     495         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    486496 
    487497      ii = ii + 1 
    … …  
    491501      ii = ii + 1 
    492502      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl),         & 
     503         &      v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) ,                                                         & 
    493504         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 
    494505         &      STAT=ierr(ii) ) 
    … …  
    505516      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      &  
    506517         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      & 
    507          &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), & 
     518         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 
    508519         &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 
    509520 

  • NEMO/trunk/src/ICE/icectl.F90

    r13601 r14005  
    8585      !! 
    8686      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
    87          &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 
     87         &          zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 
     88         &          zdiag_eimin, zdiag_esmin, zdiag_simin 
    8889      REAL(wp) ::   zvtrp, zetrp 
    8990      REAL(wp) ::   zarea 
    … …  
    9293      IF( icount == 0 ) THEN 
    9394 
    94          pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 
     95         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 
    9596         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
    9697         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
    … …  
    112113 
    113114         ! -- mass diag -- ! 
    114          zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_Dt_ice       & 
     115         zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t )      & 
     116            &            - pdiag_v ) * r1_Dt_ice                                                                          & 
    115117            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
    116118            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
    … …  
    132134 
    133135         ! -- min/max diag -- ! 
    134          zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    135          zdiag_vmin  = glob_min( 'icectl', v_i ) 
    136          zdiag_amin  = glob_min( 'icectl', a_i ) 
    137          zdiag_smin  = glob_min( 'icectl', sv_i ) 
     136         zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     137         zdiag_vimin = glob_min( 'icectl', v_i  ) 
     138         zdiag_vsmin = glob_min( 'icectl', v_s  ) 
     139         zdiag_vpmin = glob_min( 'icectl', v_ip ) 
     140         zdiag_vlmin = glob_min( 'icectl', v_il ) 
     141         zdiag_aimin = glob_min( 'icectl', a_i  ) 
     142         zdiag_simin = glob_min( 'icectl', sv_i ) 
    138143         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
    139144         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    … …  
    155160               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
    156161            ! check negative values 
    157             IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin 
    158             IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin 
    159             IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin 
    160             IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin 
    161             IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin 
     162            IF( zdiag_vimin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i  < 0        = ',zdiag_vimin 
     163            IF( zdiag_vsmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_s  < 0        = ',zdiag_vsmin 
     164            IF( zdiag_vpmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_ip < 0        = ',zdiag_vpmin 
     165            IF( zdiag_vlmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_il < 0        = ',zdiag_vlmin 
     166            IF( zdiag_aimin < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i  < 0        = ',zdiag_aimin 
     167            IF( zdiag_simin < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i  < 0        = ',zdiag_simin 
     168            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i  < 0        = ',zdiag_eimin 
     169            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s  < 0        = ',zdiag_esmin 
    162170            ! check maximum ice concentration 
    163             IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
    164                &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
     171            IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     172               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_aimax 
    165173            ! check if advection scheme is conservative 
    166174            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    167                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
     175               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 
    168176            IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    169                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rdt_ice 
     177               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rDt_ice 
    170178         ENDIF 
    171179         ! 
    … …  
    193201      ! water flux 
    194202      ! -- mass diag -- ! 
    195       zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
    196          &                              + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 
     203      zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr   + wfx_sub + wfx_pnd & 
     204         &                              + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 
    197205 
    198206      ! -- salt diag -- ! 
    … …  
    200208 
    201209      ! -- heat diag -- ! 
    202       zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     210      zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
    203211      ! equivalent to this: 
    204212      !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
    … …  
    245253      IF( icount == 0 ) THEN 
    246254 
    247          pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 ) 
    248          pdiag_s = SUM( sv_i * rhoi             , dim=3 ) 
     255         pdiag_v = SUM( v_i  * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 
     256         pdiag_s = SUM( sv_i * rhoi , dim=3 ) 
    249257         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 
    250258 
    … …  
    261269 
    262270         ! -- mass diag -- ! 
    263          zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_Dt_ice                             & 
     271         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice    & 
    264272            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
    265273            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           & 
    … …  
    352360      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    !  
    353361      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    !  
     362      ! mean state 
     363      CALL iom_rstput( 0, 0, inum, 'icecon'    , SUM(a_i ,dim=3) , ktype = jp_r8 )    ! 
     364      CALL iom_rstput( 0, 0, inum, 'icevol'    , SUM(v_i ,dim=3) , ktype = jp_r8 )    ! 
     365      CALL iom_rstput( 0, 0, inum, 'snwvol'    , SUM(v_s ,dim=3) , ktype = jp_r8 )    ! 
     366      CALL iom_rstput( 0, 0, inum, 'pndvol'    , SUM(v_ip,dim=3) , ktype = jp_r8 )    ! 
     367      CALL iom_rstput( 0, 0, inum, 'lidvol'    , SUM(v_il,dim=3) , ktype = jp_r8 )    ! 
    354368       
    355369      CALL iom_close( inum ) 
    … …  
    776790      ! -- mass diag -- ! 
    777791      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
    778          &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 
     792         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rDt_ice 
    779793      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 
    780794 
    781795      ! -- salt diag -- ! 
    782       zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 
     796      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 
    783797      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 
    784798 

  • NEMO/trunk/src/ICE/icedyn.F90

    r13472 r14005  
    2929   USE lbclnk         ! lateral boundary conditions (or mpp links) 
    3030   USE timing         ! Timing 
     31   USE fldread        ! read input fields 
    3132 
    3233   IMPLICIT NONE 
    … …  
    5152   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D) 
    5253    
     54   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_icbmsk   ! structure of input grounded icebergs mask (file informations, fields read) 
     55 
    5356   !! * Substitutions 
    5457#  include "do_loop_substitute.h90" 
    … …  
    8184      ! 
    8285      ! controls 
    83       IF( ln_timing )   CALL timing_start('icedyn') 
     86      IF( ln_timing )   CALL timing_start('ice_dyn') 
    8487      ! 
    8588      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    106109      END WHERE 
    107110      ! 
     111      IF( ln_landfast_L16 ) THEN 
     112         CALL fld_read( kt, 1, sf_icbmsk ) 
     113         icb_mask(:,:) = sf_icbmsk(1)%fnow(:,:,1) 
     114      ENDIF 
    108115      ! 
    109116      SELECT CASE( nice_dyn )          !-- Set which dynamics is running 
    … …  
    172179      ! 
    173180      ! controls 
    174       IF( ln_timing )   CALL timing_stop ('icedyn') 
     181      IF( ln_timing )   CALL timing_stop ('ice_dyn') 
    175182      ! 
    176183   END SUBROUTINE ice_dyn 
    … …  
    216223      !! ** input   :   Namelist namdyn 
    217224      !!------------------------------------------------------------------- 
    218       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
     225      INTEGER ::   ios, ioptio, ierror   ! Local integer output status for namelist read 
     226      ! 
     227      CHARACTER(len=256) ::   cn_dir     ! Root directory for location of ice files 
     228      TYPE(FLD_N)        ::   sn_icbmsk  ! informations about the grounded icebergs field to be read 
    219229      !! 
    220230      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
    221231         &             rn_ishlat ,                                                           & 
    222          &             ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile 
     232         &             ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile, & 
     233         &             sn_icbmsk, cn_dir 
    223234      !!------------------------------------------------------------------- 
    224235      ! 
    … …  
    269280      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp 
    270281      ! 
     282      !                                      !--- allocate and fill structure for grounded icebergs mask 
     283      IF( ln_landfast_L16 ) THEN 
     284         ALLOCATE( sf_icbmsk(1), STAT=ierror ) 
     285         IF( ierror > 0 ) THEN 
     286            CALL ctl_stop( 'ice_dyn_init: unable to allocate sf_icbmsk structure' ) ; RETURN 
     287         ENDIF 
     288         ! 
     289         CALL fld_fill( sf_icbmsk, (/ sn_icbmsk /), cn_dir, 'ice_dyn_init',   & 
     290            &                                               'landfast ice is a function of read grounded icebergs', 'icedyn' ) 
     291         ! 
     292         ALLOCATE( sf_icbmsk(1)%fnow(jpi,jpj,1) ) 
     293         IF( sf_icbmsk(1)%ln_tint )   ALLOCATE( sf_icbmsk(1)%fdta(jpi,jpj,1,2) ) 
     294         IF( TRIM(sf_icbmsk(1)%clrootname) == 'NOT USED' ) sf_icbmsk(1)%fnow(:,:,1) = 0._wp   ! not used field  (set to 0)          
     295      ELSE 
     296         icb_mask(:,:) = 0._wp 
     297      ENDIF 
     298      !                                      !--- other init  
    271299      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 
    272300      CALL ice_dyn_rhg_init             ! set ice rheology parameters 

  • NEMO/trunk/src/ICE/icedyn_adv_pra.F90

    r13970 r14005  
    156156 
    157157         ! diagnostics 
    158          zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     158         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     159            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
    159160         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
    160161         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
    … …  
    178179               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    179180            END DO 
    180             IF ( ln_pnd_LEV ) THEN 
     181            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    181182               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
    182183               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
    … …  
    214215            END DO 
    215216            ! 
    216             IF ( ln_pnd_LEV ) THEN 
     217            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    217218               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    218219               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    … …  
    249250                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    250251            END DO 
    251             IF ( ln_pnd_LEV ) THEN 
     252            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    252253               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    253254               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    … …  
    278279         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
    279280            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
    280          IF ( ln_pnd_LEV ) THEN 
     281         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    281282            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
    282283               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
    … …  
    302303               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    303304            END DO 
    304             IF ( ln_pnd_LEV ) THEN 
     305            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    305306               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    306307               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    … …  
    320321         ! 
    321322         ! --- diagnostics --- ! 
    322          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     323         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     324            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    323325            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    324326         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
    … …  
    769771               !                               ! -- check h_ip -- ! 
    770772               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    771                IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     773               IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    772774                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    773775                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    … …  
    10151017            END DO 
    10161018            ! 
    1017             IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     1019            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction 
    10181020               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
    10191021                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 
    … …  
    10571059            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    10581060            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    1059             IF( ln_pnd_LEV ) THEN 
     1061            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    10601062               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
    10611063               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
    … …  
    11351137         END DO 
    11361138         ! 
    1137          IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
     1139         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction 
    11381140            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    11391141            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 

  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r13633 r14005  
    182182 
    183183         ! diagnostics 
    184          zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     184         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     185            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
    185186         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
    186187         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
    … …  
    338339         ! 
    339340         !== melt ponds ==! 
    340          IF ( ln_pnd_LEV ) THEN 
     341         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    341342            ! concentration 
    342343            zamsk = 1._wp 
    … …  
    358359 
    359360         ! --- Lateral boundary conditions --- ! 
    360          IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     361         IF    ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 
    361362            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    362363               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
    363          ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     364         ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 
    364365            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    365366               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
    … …  
    379380         ! 
    380381         ! --- diagnostics --- ! 
    381          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     382         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     383            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    382384            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    383385         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
    … …  
    14971499               !                               ! -- check h_ip -- ! 
    14981500               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1499                IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1501               IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    15001502                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    15011503                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 

  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r13999 r14005  
    578578               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft  
    579579 
    580                IF ( ln_pnd_LEV ) THEN 
     580               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    581581                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    582582                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    … …  
    615615               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji) 
    616616               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    617                IF ( ln_pnd_LEV ) THEN 
     617               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    618618                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    619619                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    … …  
    712712                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
    713713                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    714                   IF ( ln_pnd_LEV ) THEN 
     714                  IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    715715                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    716716                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 

  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r13982 r14005  
    326326            zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    327327            ! ice-bottom stress at U points 
    328             zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 
     328            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    329329            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    330330            ! ice-bottom stress at V points 
    331             zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 
     331            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    332332            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    333333            ! ice_bottom stress at T points 
    334             zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 
     334            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0 
    335335            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    336336         END_2D 

  • NEMO/trunk/src/ICE/iceistate.F90

    r13472 r14005  
    426426      ! 4) Snow-ice mass (case ice is fully embedded) 
    427427      !---------------------------------------------- 
    428       snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     428      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3  )   ! snow+ice mass 
    429429      snwice_mass_b(:,:) = snwice_mass(:,:) 
    430430      ! 

  • NEMO/trunk/src/ICE/iceitd.F90

    r13618 r14005  
    2929   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    3030   USE prtctl         ! Print control 
     31   USE timing         ! Timing 
    3132 
    3233   IMPLICIT NONE 
    … …  
    8788      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories 
    8889      !!------------------------------------------------------------------ 
     90      IF( ln_timing )   CALL timing_start('iceitd_rem') 
    8991 
    9092      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'  
    … …  
    315317            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    316318               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    317                IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
     319               IF( ln_pnd_LEV .OR. ln_pnd_TOPO )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    318320               h_i_1d(ji) = rn_himin 
    319321            ENDIF 
    … …  
    328330      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    329331      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     332      IF( ln_timing    )   CALL timing_stop ('iceitd_rem') 
    330333      ! 
    331334   END SUBROUTINE ice_itd_rem 
    … …  
    486489               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    487490               !   
    488                IF ( ln_pnd_LEV ) THEN 
     491               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    489492                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction 
    490493                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 
    491494                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    492495                  !                                               
    493                   ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
     496                  ztrans          = v_ip_2d(ji,jl1) * zworkv(ji)     ! Pond volume 
    494497                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    495498                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    496499                  ! 
    497500                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume 
    498                      ztrans          = v_il_2d(ji,jl1) * zworka(ji) 
     501                     ztrans          = v_il_2d(ji,jl1) * zworkv(ji) 
    499502                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 
    500503                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans 
    … …  
    606609      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    607610      !!------------------------------------------------------------------ 
     611      IF( ln_timing )   CALL timing_start('iceitd_reb') 
    608612      ! 
    609613      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
    … …  
    635639               jdonor(ji,jl)  = jl  
    636640               ! how much of a_i you send in cat sup is somewhat arbitrary 
    637                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    638                !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    639                !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    640                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    641                !!          zdaice(ji,jl)  = a_i_1d(ji) 
    642                !!          zdvice(ji,jl)  = v_i_1d(ji) 
    643                !!clem: these are from UCL and work ok 
    644                zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    645                zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     641               ! these are from CICE => transfer everything 
     642               !!zdaice(ji,jl)  = a_i_1d(ji) 
     643               !!zdvice(ji,jl)  = v_i_1d(ji) 
     644               ! these are from LLN => transfer only half of the category 
     645               zdaice(ji,jl)  =                       0.5_wp  * a_i_1d(ji) 
     646               zdvice(ji,jl)  = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 
    646647            END DO 
    647648            ! 
    … …  
    686687      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    687688      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     689      IF( ln_timing    )   CALL timing_stop ('iceitd_reb') 
    688690      ! 
    689691   END SUBROUTINE ice_itd_reb 

  • NEMO/trunk/src/ICE/icesbc.F90

    r13472 r14005  
    6262      !!------------------------------------------------------------------- 
    6363      ! 
    64       IF( ln_timing )   CALL timing_start('ice_sbc') 
     64      IF( ln_timing )   CALL timing_start('icesbc') 
    6565      ! 
    6666      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    8989      ENDIF 
    9090      ! 
    91       IF( ln_timing )   CALL timing_stop('ice_sbc') 
     91      IF( ln_timing )   CALL timing_stop('icesbc') 
    9292      ! 
    9393   END SUBROUTINE ice_sbc_tau 
    … …  
    122122      !!-------------------------------------------------------------------- 
    123123      ! 
    124       IF( ln_timing )   CALL timing_start('ice_sbc_flx') 
     124      IF( ln_timing )   CALL timing_start('icesbc') 
    125125 
    126126      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    176176      ENDIF 
    177177      ! 
    178       IF( ln_timing )   CALL timing_stop('ice_sbc_flx') 
     178      IF( ln_timing )   CALL timing_stop('icesbc') 
    179179      ! 
    180180   END SUBROUTINE ice_sbc_flx 

  • NEMO/trunk/src/ICE/icestp.F90

    r13970 r14005  
    121121      !!---------------------------------------------------------------------- 
    122122      ! 
    123       IF( ln_timing )   CALL timing_start('ice_stp') 
     123      IF( ln_timing )   CALL timing_start('icestp') 
    124124      ! 
    125125      !                                      !-----------------------! 
    … …  
    215215!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    216216      ! 
    217       IF( ln_timing )   CALL timing_stop('ice_stp') 
     217      IF( ln_timing )   CALL timing_stop('icestp') 
    218218      ! 
    219219   END SUBROUTINE ice_stp 
    … …  
    373373      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume 
    374374      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume 
     375      v_ip_b(:,:,:)   = v_ip(:,:,:)     ! pond volume 
     376      v_il_b(:,:,:)   = v_il(:,:,:)     ! pond lid volume 
    375377      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content 
    376378      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy 
    … …  
    432434         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
    433435         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
     436         diag_aice(ji,jj) = 0._wp ;   diag_vpnd(ji,jj) = 0._wp 
    434437 
    435438         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    … …  
    457460            qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T) 
    458461            qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs 
     462            ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux) 
     463            dh_i_sum_2d(ji,jj,jl) = 0._wp 
     464            dh_s_mlt_2d(ji,jj,jl) = 0._wp 
    459465         END_2D 
    460466      ENDDO 
    … …  
    485491         diag_vsnw(:,:) = diag_vsnw(:,:) & 
    486492            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
     493         diag_vpnd(:,:) = diag_vpnd(:,:) & 
     494            &             + SUM(     v_ip + v_il          - v_ip_b - v_il_b                , dim=3 ) * r1_Dt_ice * rhow 
    487495         ! 
    488496         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 

  • NEMO/trunk/src/ICE/icethd.F90

    r13643 r14005  
    166166         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
    167167 
     168         ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 
     169         ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 
     170         IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 
     171            zqfr               = 0._wp 
     172            zqfr_pos           = 0._wp 
     173            qsb_ice_bot(ji,jj) = 0._wp 
     174         ENDIF 
     175         ! 
    168176         ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
    169177         !     qlead is the energy received from the atm. in the leads. 
    … …  
    239247            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
    240248                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    241                               CALL ice_thd_pnd                          ! Melt ponds formation 
    242249                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    243250            ENDIF 
    … …  
    260267      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    261268      !                    
     269      IF ( ln_pnd .AND. ln_icedH ) & 
     270         &                    CALL ice_thd_pnd                      ! --- Melt ponds  
     271      ! 
    262272      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
    263273      ! 
    … …  
    266276                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
    267277      ! 
    268       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation      
     278      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! ice natural aging incrementation      
    269279      ! 
    270280      ! convergence tests 
    … …  
    377387            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
    378388         END DO 
    379          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
    380          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    381          CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) ) 
    382389         ! 
    383390         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
    … …  
    409416         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          ) 
    410417         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          ) 
    411          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          ) 
    412418         ! 
    413419         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          ) 
    … …  
    464470         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 
    465471         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 
    466          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    467          v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
    468472         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 
    469473          
    … …  
    483487            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
    484488         END DO 
    485          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
    486          CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    487          CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) ) 
    488489         ! 
    489490         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
    … …  
    501502         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        ) 
    502503         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        ) 
    503          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    504504         ! 
    505505         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        ) 
    … …  
    529529         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) ) 
    530530         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) ) 
     531         ! Melt ponds 
     532         CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum  (1:npti) , dh_i_sum_2d(:,:,kl) ) 
     533         CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt  (1:npti) , dh_s_mlt_2d(:,:,kl) ) 
    531534         ! SIMIP diagnostics          
    532535         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) ) 
    … …  
    537540         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 
    538541         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 
    539          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 
    540          CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) ) 
    541542         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 
    542543         ! check convergence of heat diffusion scheme 

  • NEMO/trunk/src/ICE/icethd_dh.F90

    r13643 r14005  
    5555      !!                - Snow ice formation 
    5656      !! 
     57      !! ** Note     :  h=max(0,h+dh) are often used to ensure positivity of h. 
     58      !!                very small negative values can occur otherwise (e.g. -1.e-20) 
     59      !! 
    5760      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 
    5861      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646    
    … …  
    7982      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
    8083 
    81       REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
    8284      REAL(wp), DIMENSION(jpij) ::   zq_top      ! heat for surface ablation                   (J.m-2) 
    8385      REAL(wp), DIMENSION(jpij) ::   zq_bot      ! heat for bottom ablation                    (J.m-2) 
    … …  
    8587      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    8688      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2) 
    87  
    88       REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt  
    89       REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation  
    90       REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation 
    91  
    92       REAL(wp), DIMENSION(jpij,nlay_s) ::   zh_s      ! snw layer thickness 
    93       REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    94       REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
    95       INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    96  
     89      REAL(wp), DIMENSION(jpij) ::   zdeltah      
    9790      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing 
     91 
     92      INTEGER , DIMENSION(jpij,nlay_i)     ::   icount    ! number of layers vanishing by melting  
     93      REAL(wp), DIMENSION(jpij,0:nlay_i+1) ::   zh_i      ! ice layer thickness (m) 
     94      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   zh_s      ! snw layer thickness (m) 
     95      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   ze_s      ! snw layer enthalpy (J.m-3) 
    9896 
    9997      REAL(wp) ::   zswitch_sal 
    … …  
    108106      END SELECT 
    109107 
    110       ! initialize layer thicknesses and enthalpies 
     108      ! initialize ice layer thicknesses and enthalpies 
     109      eh_i_old(1:npti,0:nlay_i+1) = 0._wp 
    111110      h_i_old (1:npti,0:nlay_i+1) = 0._wp 
    112       eh_i_old(1:npti,0:nlay_i+1) = 0._wp 
     111      zh_i    (1:npti,0:nlay_i+1) = 0._wp 
    113112      DO jk = 1, nlay_i 
    114113         DO ji = 1, npti 
     114            eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 
    115115            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 
    116             eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 
     116            zh_i    (ji,jk) = h_i_1d(ji) * r1_nlay_i 
     117         END DO 
     118      END DO 
     119      ! 
     120      ! initialize snw layer thicknesses and enthalpies 
     121      zh_s(1:npti,0) = 0._wp 
     122      ze_s(1:npti,0) = 0._wp 
     123      DO jk = 1, nlay_s 
     124         DO ji = 1, npti 
     125            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
     126            ze_s(ji,jk) = e_s_1d(ji,jk) 
    117127         END DO 
    118128      END DO 
    … …  
    141151         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji)  
    142152         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 
    143       END DO 
    144  
    145       ! Ice and snow layer thicknesses 
    146       !------------------------------- 
    147       DO jk = 1, nlay_i 
    148          DO ji = 1, npti 
    149             zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i 
    150          END DO 
    151       END DO 
    152  
    153       DO jk = 1, nlay_s 
    154          DO ji = 1, npti 
    155             zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
    156          END DO 
    157153      END DO 
    158154 
    … …  
    167163         DO ji = 1, npti 
    168164            IF( t_s_1d(ji,jk) > rt0 ) THEN 
    169                hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0 
    170                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux 
     165               hfx_res_1d    (ji) = hfx_res_1d    (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0 
     166               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux 
    171167               ! updates 
    172                dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
    173                h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk) 
     168               dh_s_mlt(ji)    =             dh_s_mlt(ji) - zh_s(ji,jk) 
     169               h_s_1d  (ji)    = MAX( 0._wp, h_s_1d  (ji) - zh_s(ji,jk) ) 
    174170               zh_s    (ji,jk) = 0._wp 
    175                e_s_1d  (ji,jk) = 0._wp 
    176                t_s_1d  (ji,jk) = rt0 
     171               ze_s    (ji,jk) = 0._wp 
    177172            END IF 
    178173         END DO 
    … …  
    181176      ! Snow precipitation 
    182177      !------------------- 
    183       CALL ice_var_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
    184  
    185       zdeltah(1:npti,:) = 0._wp 
     178      CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
     179 
    186180      DO ji = 1, npti 
    187181         IF( sprecip_1d(ji) > 0._wp ) THEN 
     182            zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji)   ! thickness of precip 
     183            ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) )                              ! enthalpy of the precip (>0, J.m-3) 
    188184            ! 
    189             ! --- precipitation --- 
    190             zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji)   ! thickness change 
    191             zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
     185            hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2) 
     186            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos       * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice   ! mass flux, <0 
    192187            ! 
    193             hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2) 
    194             wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0 
    195              
    196             ! --- melt of falling snow --- 
    197             rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
    198             zdeltah       (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change 
    199             zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting  
    200             hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat used to melt snow (W.m-2, >0) 
    201             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0 
    202              
    203             ! updates available heat + precipitations after melting 
    204             dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1) 
    205             zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )       
    206             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    207              
    208188            ! update thickness 
    209             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 
    210             ! 
    211          ELSE 
    212             ! 
    213             zdh_s_pre(ji) = 0._wp 
    214             zqprec   (ji) = 0._wp 
    215             ! 
     189            h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 
    216190         ENDIF 
    217       END DO 
    218  
    219       ! recalculate snow layers 
    220       DO jk = 1, nlay_s 
    221          DO ji = 1, npti 
    222             zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
    223          END DO 
    224191      END DO 
    225192 
    226193      ! Snow melting 
    227194      ! ------------ 
    228       ! If heat still available (zq_top > 0), then melt more snow 
    229       zdeltah(1:npti,:) = 0._wp 
    230       zdh_s_mel(1:npti) = 0._wp 
    231       DO jk = 1, nlay_s 
     195      ! If heat still available (zq_top > 0) 
     196      ! then all snw precip has been melted and we need to melt more snow 
     197      DO jk = 0, nlay_s 
    232198         DO ji = 1, npti 
    233199            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 
    234200               ! 
    235                rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 
    236                zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change 
    237                zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
    238                zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
    239                 
    240                hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0) 
    241                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip) 
     201               rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) ) 
     202               zdum    = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 )   ! thickness change 
     203               zdum    = MAX( zdum , - zh_s(ji,jk) )                           ! bound melting 
     204                
     205               hfx_snw_1d    (ji) = hfx_snw_1d    (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0) 
     206               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_Dt_ice   ! snow melting only = water into the ocean 
    242207                
    243208               ! updates available heat + thickness 
    244                dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk) 
    245                zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    246                h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    247                zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 
     209               dh_s_mlt(ji)    =              dh_s_mlt(ji)    + zdum 
     210               zq_top  (ji)    = MAX( 0._wp , zq_top  (ji)    + zdum * ze_s(ji,jk) ) 
     211               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d  (ji)    + zdum ) 
     212               zh_s    (ji,jk) = MAX( 0._wp , zh_s    (ji,jk) + zdum ) 
     213!!$               IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp 
    248214               ! 
    249215            ENDIF 
    … …  
    255221      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    256222      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 
    257       zdeltah(1:npti,:) = 0._wp 
     223      zdeltah   (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 
     224      zevap_rema(1:npti) = 0._wp 
    258225      DO ji = 1, npti 
    259226         IF( evap_ice_1d(ji) > 0._wp ) THEN 
     227            zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, < 0             
     228            zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on) 
     229         ENDIF 
     230      END DO 
     231       
     232      DO jk = 0, nlay_s 
     233         DO ji = 1, npti 
     234            zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 
    260235            ! 
    261             zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rDt_ice ) 
    262             zevap_rema(ji)   = evap_ice_1d(ji) * rDt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on) 
    263             zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    264              
    265             hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 
    266                &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
    267                &                 * a_i_1d(ji) * r1_Dt_ice 
    268             wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation 
    269              
    270             ! new snow thickness 
    271             h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 
    272             ! update precipitations after sublimation and correct sublimation 
    273             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    274             zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
    275             ! 
    276          ELSE 
    277             ! 
    278             zdh_s_sub (ji) = 0._wp 
    279             zevap_rema(ji) = 0._wp 
    280             ! 
    281          ENDIF 
    282       END DO 
    283        
    284       ! --- Update snow diags --- ! 
    285       DO ji = 1, npti 
    286          dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    287       END DO 
    288  
    289       ! Update temperature, energy 
    290       !--------------------------- 
    291       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
    292       DO jk = 1, nlay_s 
    293          DO ji = 1,npti 
    294             rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 
    295             e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
    296               &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
    297               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 
    298          END DO 
    299       END DO 
    300        
     236            hfx_sub_1d    (ji) = hfx_sub_1d    (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! Heat flux of snw that sublimates [W.m-2], < 0 
     237            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_Dt_ice  ! Mass flux by sublimation 
     238 
     239            ! update thickness 
     240            h_s_1d(ji)    = MAX( 0._wp , h_s_1d(ji)    + zdum ) 
     241            zh_s  (ji,jk) = MAX( 0._wp , zh_s  (ji,jk) + zdum ) 
     242!!$            IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp 
     243 
     244            ! update sublimation left 
     245            zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) 
     246         END DO 
     247      END DO 
     248 
     249      !       
    301250      !                       ! ============ ! 
    302251      !                       !     Ice      ! 
    … …  
    305254      ! Surface ice melting  
    306255      !-------------------- 
    307       zdeltah(1:npti,:) = 0._wp ! important 
    308256      DO jk = 1, nlay_i 
    309257         DO ji = 1, npti 
    … …  
    313261 
    314262               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]        
    315                zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    316                                                                       ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    317                zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
    318                                                                       ! this should normally not happen, but sometimes, heat diffusion leads to this 
    319                zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0 
    320                           
    321                dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting 
    322                 
    323                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    324  
    325                hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    326                !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    327                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice    ! Salt flux 
    328                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    329                wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                 ! Mass flux 
    330  
     263               zdE            =   0._wp                               ! Specific enthalpy difference (J/kg, <0) 
     264               !                                                          set up at 0 since no energy is needed to melt water...(it is already melted) 
     265               zdum           = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     266               !                                                          this should normally not happen, but sometimes, heat diffusion leads to this 
     267               zfmdt          = - zdum * rhoi                         ! Recompute mass flux [kg/m2, >0] 
     268               ! 
     269               dh_i_itm(ji)   = dh_i_itm(ji) + zdum                   ! Cumulate internal melting 
     270               ! 
     271               hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux to the ocean [W.m-2], <0 
     272               !                                                                                          ice enthalpy zEi is "sent" to the ocean 
     273               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice    ! Mass flux 
     274               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice    ! Salt flux 
     275               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok 
    331276            ELSE                                        !-- Surface melting 
    332277                
    … …  
    337282               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0] 
    338283                
    339                zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    340                 
    341                zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    342                 
    343                zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    344                 
    345                dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    346                 
    347                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
     284               zdum          = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
     285                
     286               zdum           = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     287 
     288               zq_top(ji)     = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat 
     289                
     290               dh_i_sum(ji)   = dh_i_sum(ji) + zdum                   ! Cumulate surface melt 
     291                
     292               zfmdt          = - rhoi * zdum                         ! Recompute mass flux [kg/m2, >0] 
    348293                
    349294               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    350295                
    351                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice    ! Salt flux >0 
    352                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    353                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux [W.m-2], < 0 
    354                hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], > 0   
    355                !  
    356                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                 ! Mass flux 
    357                 
     296               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux [W.m-2], < 0 
     297               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_Dt_ice    ! Heat flux used in this process [W.m-2], > 0   
     298               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice    ! Mass flux 
     299               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice    ! Salt flux >0 
     300               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok)  
    358301            END IF 
    359              
     302            ! update thickness 
     303            zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     304            h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     305            ! 
     306            ! update heat content (J.m-2) and layer thickness 
     307            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     308            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
     309            ! 
     310            ! 
    360311            ! Ice sublimation 
    361312            ! --------------- 
    362             zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 
    363             zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
    364             dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
    365              
    366             sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice  ! Salt flux >0 
    367             !                                                                                          clem: flux is sent to the ocean for simplicity 
    368             !                                                                                                but salt should remain in the ice except 
    369             !                                                                                                if all ice is melted. => must be corrected 
    370             hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice      ! Heat flux [W.m-2], < 0 
    371  
    372             wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice           ! Mass flux > 0 
    373  
    374             ! update remaining mass flux 
    375             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi 
    376              
     313            zdum               = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) 
     314            ! 
     315            hfx_sub_1d(ji)     = hfx_sub_1d(ji)     + e_i_1d(ji,jk) * zdum              * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 
     316            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi          * zdum              * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0 
     317            sfx_sub_1d(ji)     = sfx_sub_1d(ji)     - rhoi          * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0 
     318            !                                                                                                      clem: flux is sent to the ocean for simplicity 
     319            !                                                                                                            but salt should remain in the ice except 
     320            !                                                                                                            if all ice is melted. => must be corrected 
     321            ! update remaining mass flux and thickness 
     322            zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi             
     323            zh_i(ji,jk)    = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     324            h_i_1d(ji)     = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     325            dh_i_sub(ji)   = dh_i_sub(ji) + zdum 
     326 
     327            ! update heat content (J.m-2) and layer thickness 
     328            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     329            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
     330 
    377331            ! record which layers have disappeared (for bottom melting)  
    378332            !    => icount=0 : no layer has vanished 
    379333            !    => icount=5 : 5 layers have vanished 
    380             rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     334            rswitch       = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) )  
    381335            icount(ji,jk) = NINT( rswitch ) 
    382             zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    383336                         
    384             ! update heat content (J.m-2) and layer thickness 
    385             eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    386             h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    387          END DO 
    388       END DO 
    389        
    390       ! update ice thickness 
    391       DO ji = 1, npti 
    392          h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 
    393       END DO 
    394  
     337         END DO 
     338      END DO 
     339       
    395340      ! remaining "potential" evap is sent to ocean 
    396341      DO ji = 1, npti 
    … …  
    430375                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
    431376 
    432                s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
    433  
    434                ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C) 
    435  
    436                zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    437                 
    438                zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
    439                   &            - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 
    440  
    441                zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
    442  
    443                zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
    444  
    445                dh_i_bog(ji)  = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
     377               s_i_new(ji)    = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
     378 
     379               ztmelts        = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C) 
     380 
     381               zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
     382                
     383               zEi            = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     384                  &             - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 
     385 
     386               zEw            = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     387 
     388               zdE            = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
     389 
     390               dh_i_bog(ji)   = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
    446391                
    447392            END DO 
    448393            ! Contribution to Energy and Salt Fluxes                                     
    449             zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
     394            zfmdt = - rhoi * dh_i_bog(ji)                                                              ! Mass flux x time step (kg/m2, < 0) 
    450395             
    451             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], >0 
    452             hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], <0 
    453              
    454             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice     ! Salt flux, <0 
    455  
    456             wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                   ! Mass flux, <0 
     396            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt                      * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], >0 
     397            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE  * zfmdt                      * a_i_1d(ji) * r1_Dt_ice   ! Heat flux used in this process [W.m-2], <0           
     398            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji)               * a_i_1d(ji) * r1_Dt_ice   ! Mass flux, <0 
     399            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux, <0 
     400 
     401            ! update thickness 
     402            zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji) 
     403            h_i_1d(ji)        = h_i_1d(ji)        + dh_i_bog(ji) 
    457404 
    458405            ! update heat content (J.m-2) and layer thickness 
    … …  
    466413      ! Ice Basal melt 
    467414      !--------------- 
    468       zdeltah(1:npti,:) = 0._wp ! important 
    469415      DO jk = nlay_i, 1, -1 
    470416         DO ji = 1, npti 
    … …  
    475421               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    476422 
    477                   zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
    478                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    479                                                                     !    set up at 0 since no energy is needed to melt water...(it is already melted) 
    480                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
    481                                                                     ! this should normally not happen, but sometimes, heat diffusion leads to this 
    482  
    483                   dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk) 
    484  
    485                   zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0 
    486  
    487                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    488                   !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    489                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice    ! Salt flux 
    490                   !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    491                   wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                 ! Mass flux 
    492  
    493                   ! update heat content (J.m-2) and layer thickness 
    494                   eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    495                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    496  
     423                  zEi            = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
     424                  zdE            = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
     425                  !                                                  set up at 0 since no energy is needed to melt water...(it is already melted) 
     426                  zdum           = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     427                  !                                                  this should normally not happen, but sometimes, heat diffusion leads to this 
     428                  dh_i_itm (ji)  = dh_i_itm(ji) + zdum 
     429                  ! 
     430                  zfmdt          = - zdum * rhoi                 ! Mass flux x time step > 0 
     431                  ! 
     432                  hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], <0 
     433                  !                                                                                         ice enthalpy zEi is "sent" to the ocean 
     434                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice   ! Mass flux 
     435                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux 
     436                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok 
    497437               ELSE                                        !-- Basal melting 
    498438 
    499                   zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0) 
    500                   zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
    501                   zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
    502  
    503                   zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0) 
    504  
    505                   zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change 
    506  
    507                   zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
     439                  zEi            = - e_i_1d(ji,jk) * r1_rhoi                       ! Specific enthalpy of melting ice (J/kg, <0) 
     440                  zEw            = rcp * ztmelts                                   ! Specific enthalpy of meltwater (J/kg, <0) 
     441                  zdE            = zEi - zEw                                       ! Specific enthalpy difference   (J/kg, <0) 
     442 
     443                  zfmdt          = - zq_bot(ji) / zdE                              ! Mass flux x time step (kg/m2, >0) 
     444 
     445                  zdum           = - zfmdt * r1_rhoi                               ! Gross thickness change 
     446 
     447                  zdum           = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) )       ! bound thickness change 
    508448                   
    509                   zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors 
    510  
    511                   dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt 
    512  
    513                   zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0 
    514  
    515                   zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
    516  
    517                   hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0   
    518                   hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat used in this process [W.m-2], >0   
    519  
    520                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
    521                   !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    522                    
    523                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                 ! Mass flux 
    524  
    525                   ! update heat content (J.m-2) and layer thickness 
    526                   eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    527                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
     449                  zq_bot(ji)     = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors 
     450 
     451                  dh_i_bom(ji)   = dh_i_bom(ji) + zdum                             ! Update basal melt 
     452 
     453                  zfmdt          = - zdum * rhoi                                   ! Mass flux x time step > 0 
     454 
     455                  zQm            = zfmdt * zEw                                     ! Heat exchanged with ocean 
     456 
     457                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat flux to the ocean [W.m-2], <0   
     458                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_Dt_ice   ! Heat used in this process [W.m-2], >0 
     459                  wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_Dt_ice   ! Mass flux 
     460                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice   ! Salt flux 
     461                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok 
    528462               ENDIF 
    529             
     463               ! update thickness 
     464               zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     465               h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     466               ! 
     467               ! update heat content (J.m-2) and layer thickness 
     468               eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     469               h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
    530470            ENDIF 
    531471         END DO 
    532472      END DO 
    533473 
    534       ! Update temperature, energy 
    535       ! -------------------------- 
    536       DO ji = 1, npti 
    537          h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 
    538       END DO   
    539  
    540       ! If heat still available then melt more snow 
    541       !------------------------------------------- 
    542       zdeltah(1:npti,:) = 0._wp ! important 
    543       DO ji = 1, npti 
    544          zq_rema (ji)   = zq_top(ji) + zq_bot(ji)  
    545          rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow 
    546          rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 
    547          zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 
    548          zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
    549          dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1) 
    550          h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1) 
    551          
    552          zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
    553          hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice   ! Heat used to melt snow, W.m-2 (>0) 
    554          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice       ! Mass flux 
    555          dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    556          !     
    557          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    558          !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 
    559  
    560          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
    561       END DO 
    562  
    563       ! 
     474      ! Remove snow if ice has melted entirely 
     475      ! -------------------------------------- 
     476      DO jk = 0, nlay_s 
     477         DO ji = 1,npti 
     478            IF( h_i_1d(ji) == 0._wp ) THEN 
     479               ! mass & energy loss to the ocean 
     480               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice  ! heat flux to the ocean [W.m-2], < 0 
     481               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice  ! mass flux 
     482 
     483               ! update thickness and energy 
     484               h_s_1d(ji)    = 0._wp 
     485               ze_s  (ji,jk) = 0._wp 
     486               zh_s  (ji,jk) = 0._wp 
     487            ENDIF 
     488         END DO 
     489      END DO 
     490       
     491      ! Snow load on ice 
     492      ! ----------------- 
     493      ! When snow load exceeds Archimede's limit and sst is positive, 
     494      ! snow-ice formation (next bloc) can lead to negative ice enthalpy. 
     495      ! Therefore we consider here that this excess of snow falls into the ocean 
     496      zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos 
     497      DO jk = 0, nlay_s 
     498         DO ji = 1, npti 
     499            IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN 
     500               ! snow layer thickness that falls into the ocean 
     501               zdum = MIN( zdeltah(ji) , zh_s(ji,jk) ) 
     502               ! mass & energy loss to the ocean 
     503               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! heat flux to the ocean [W.m-2], < 0 
     504               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zdum * a_i_1d(ji) * r1_Dt_ice  ! mass flux 
     505               ! update thickness and energy 
     506               h_s_1d(ji)    = MAX( 0._wp, h_s_1d(ji)  - zdum ) 
     507               zh_s  (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) 
     508               ! update snow thickness that still has to fall 
     509               zdeltah(ji)   = MAX( 0._wp, zdeltah(ji) - zdum ) 
     510            ENDIF 
     511         END DO 
     512      END DO 
     513       
    564514      ! Snow-Ice formation 
    565515      ! ------------------ 
    566       ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    567       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
     516      ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,  
     517      ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice) 
    568518      z1_rho = 1._wp / ( rhos+rho0-rhoi ) 
     519      zdeltah(1:npti) = 0._wp 
    569520      DO ji = 1, npti 
    570521         ! 
    571          dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 
     522         dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 
    572523 
    573524         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji) 
    … …  
    579530         zQm            = zfmdt * zEw  
    580531          
    581          hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 
    582  
    583          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 
     532         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw        * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux 
     533         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux 
    584534 
    585535         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
    586536         IF( nn_icesal /= 2 )  THEN 
    587             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_Dt_ice  & ! put back sss_m     into the ocean 
    588                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice     ! and get  rn_icesal from the ocean  
     537            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt                 * a_i_1d(ji) * r1_Dt_ice  & ! put back sss_m     into the ocean 
     538               &                            - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice     ! and get  rn_icesal from the ocean  
    589539         ENDIF 
    590540 
    591541         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    592          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 
    593          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 
     542         wfx_sni_1d    (ji) = wfx_sni_1d    (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice 
     543         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice 
     544 
     545         ! update thickness 
     546         zh_i(ji,0)  = zh_i(ji,0) + dh_snowice(ji) 
     547         zdeltah(ji) =              dh_snowice(ji) 
    594548 
    595549         ! update heat content (J.m-2) and layer thickness 
    596          eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw 
    597550         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    598           
    599       END DO 
    600  
    601       ! 
    602       ! Update temperature, energy 
    603       ! -------------------------- 
    604       DO ji = 1, npti 
    605          rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )  
    606          t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 
    607       END DO 
    608  
     551         eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw           ! 1st part (sea water enthalpy) 
     552 
     553      END DO 
     554      ! 
     555      DO jk = nlay_s, 0, -1   ! flooding of snow starts from the base 
     556         DO ji = 1, npti 
     557            zdum           = MIN( zdeltah(ji), zh_s(ji,jk) )     ! amount of snw that floods, > 0 
     558            zh_s(ji,jk)    = MAX( 0._wp, zh_s(ji,jk) - zdum )    ! remove some snow thickness 
     559            eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy) 
     560            ! update dh_snowice 
     561            zdeltah(ji)    = MAX( 0._wp, zdeltah(ji) - zdum ) 
     562         END DO 
     563      END DO 
     564      ! 
     565      ! 
     566!!$      ! --- Update snow diags --- ! 
     567!!$      !!clem: this is wrong. dh_s_tot is not used anyway 
     568!!$      DO ji = 1, npti 
     569!!$         dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji) 
     570!!$      END DO 
     571      ! 
     572      ! 
     573      ! Remapping of snw enthalpy on a regular grid 
     574      !-------------------------------------------- 
     575      CALL snw_ent( zh_s, ze_s, e_s_1d ) 
     576       
     577      ! recalculate t_s_1d from e_s_1d 
    609578      DO jk = 1, nlay_s 
    610579         DO ji = 1,npti 
    611             ! where there is no ice or no snow 
    612             rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) ) 
    613             ! mass & energy loss to the ocean 
    614             hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 
    615                &                              ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice )  ! heat flux to the ocean [W.m-2], < 0 
    616             wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 
    617                &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice )  ! mass flux 
    618             ! update energy (mass is updated in the next loop) 
    619             e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    620             ! recalculate t_s_1d from e_s_1d 
    621             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 
    622          END DO 
    623       END DO 
     580            IF( h_s_1d(ji) > 0._wp ) THEN 
     581               t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 
     582            ELSE 
     583               t_s_1d(ji,jk) = rt0 
     584            ENDIF 
     585         END DO 
     586      END DO 
     587 
     588      ! Note: remapping of ice enthalpy is done in icethd.F90 
    624589 
    625590      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 
    626591      WHERE( h_i_1d(1:npti) == 0._wp )    
    627          a_i_1d(1:npti) = 0._wp 
    628          h_s_1d(1:npti) = 0._wp 
     592         a_i_1d (1:npti) = 0._wp 
     593         h_s_1d (1:npti) = 0._wp 
     594         t_su_1d(1:npti) = rt0 
    629595      END WHERE 
    630       ! 
     596       
    631597   END SUBROUTINE ice_thd_dh 
    632598 
     599   SUBROUTINE snw_ent( ph_old, pe_old, pe_new ) 
     600      !!------------------------------------------------------------------- 
     601      !!               ***   ROUTINE snw_ent  *** 
     602      !! 
     603      !! ** Purpose : 
     604      !!           This routine computes new vertical grids in the snow,  
     605      !!           and consistently redistributes temperatures.  
     606      !!           Redistribution is made so as to ensure to energy conservation 
     607      !! 
     608      !! 
     609      !! ** Method  : linear conservative remapping 
     610      !!            
     611      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 
     612      !!            2) linear remapping on the new layers 
     613      !! 
     614      !! ------------ cum0(0)                        ------------- cum1(0) 
     615      !!                                    NEW      ------------- 
     616      !! ------------ cum0(1)               ==>      ------------- 
     617      !!     ...                                     ------------- 
     618      !! ------------                                ------------- 
     619      !! ------------ cum0(nlay_s+1)                 ------------- cum1(nlay_s) 
     620      !! 
     621      !! 
     622      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 
     623      !!------------------------------------------------------------------- 
     624      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   ph_old             ! old thicknesses (m) 
     625      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   pe_old             ! old enthlapies (J.m-3) 
     626      REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) ::   pe_new             ! new enthlapies (J.m-3, remapped) 
     627      ! 
     628      INTEGER  :: ji         !  dummy loop indices 
     629      INTEGER  :: jk0, jk1   !  old/new layer indices 
     630      ! 
     631      REAL(wp), DIMENSION(jpij,0:nlay_s+1) ::   zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces 
     632      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces 
     633      REAL(wp), DIMENSION(jpij)            ::   zhnew               ! new layers thicknesses 
     634      !!------------------------------------------------------------------- 
     635 
     636      !-------------------------------------------------------------------------- 
     637      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces 
     638      !-------------------------------------------------------------------------- 
     639      zeh_cum0(1:npti,0) = 0._wp  
     640      zh_cum0 (1:npti,0) = 0._wp 
     641      DO jk0 = 1, nlay_s+1 
     642         DO ji = 1, npti 
     643            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1) 
     644            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1) 
     645         END DO 
     646      END DO 
     647 
     648      !------------------------------------ 
     649      !  2) Interpolation on the new layers 
     650      !------------------------------------ 
     651      ! new layer thickesses 
     652      DO ji = 1, npti 
     653         zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s   
     654      END DO 
     655 
     656      ! new layers interfaces 
     657      zh_cum1(1:npti,0) = 0._wp 
     658      DO jk1 = 1, nlay_s 
     659         DO ji = 1, npti 
     660            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 
     661         END DO 
     662      END DO 
     663 
     664      zeh_cum1(1:npti,0:nlay_s) = 0._wp  
     665      ! new cumulative q*h => linear interpolation 
     666      DO jk0 = 1, nlay_s+1 
     667         DO jk1 = 1, nlay_s-1 
     668            DO ji = 1, npti 
     669               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 
     670                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  & 
     671                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  & 
     672                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 
     673               ENDIF 
     674            END DO 
     675         END DO 
     676      END DO 
     677      ! to ensure that total heat content is strictly conserved, set: 
     678      zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1)  
     679 
     680      ! new enthalpies 
     681      DO jk1 = 1, nlay_s 
     682         DO ji = 1, npti 
     683            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )  
     684            pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 
     685         END DO 
     686      END DO 
     687       
     688   END SUBROUTINE snw_ent 
     689 
     690    
    633691#else 
    634692   !!---------------------------------------------------------------------- 

  • NEMO/trunk/src/ICE/icethd_pnd.F90

    r13472 r14005  
    2020   USE ice1D          ! sea-ice: thermodynamics variables 
    2121   USE icetab         ! sea-ice: 1D <==> 2D transformation 
     22   USE sbc_ice        ! surface energy budget 
    2223   ! 
    2324   USE in_out_manager ! I/O manager 
     25   USE iom            ! I/O manager library 
    2426   USE lib_mpp        ! MPP library 
    2527   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    … …  
    3436   INTEGER ::              nice_pnd    ! choice of the type of pond scheme 
    3537   !                                   ! associated indices: 
    36    INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme 
    37    INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant ice pond scheme 
    38    INTEGER, PARAMETER ::   np_pndLEV = 2   ! Level ice pond scheme 
    39  
     38   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme 
     39   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme 
     40   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme 
     41   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme 
     42 
     43   !--------------------------------------------------------------------------  
     44   ! Diagnostics for pond volume per area 
     45   ! 
     46   ! dV/dt = mlt + drn + lid + rnf 
     47   ! mlt   = input from surface melting 
     48   ! drn   = drainage through brine network 
     49   ! lid   = lid growth & melt 
     50   ! rnf   = runoff (water directly removed out of surface melting + overflow) 
     51   ! 
     52   ! In topo mode, the pond water lost because it is in the snow is not included in the budget 
     53   ! In level mode, all terms are incorporated 
     54   ! 
     55   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_mlt       ! meltwater pond volume input      [kg/m2/s] 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_drn       ! pond volume lost by drainage     [-] 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_lid       ! exchange with lid / refreezing   [-] 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_rnf       ! meltwater pond lost to runoff    [-]       
     59   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_mlt_1d    ! meltwater pond volume input      [-] 
     60   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_drn_1d    ! pond volume lost by drainage     [-] 
     61   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_lid_1d    ! exchange with lid / refreezing   [-] 
     62   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_rnf_1d    ! meltwater pond lost to runoff    [-] 
     63 
     64   !! * Substitutions 
     65#  include "do_loop_substitute.h90" 
    4066   !!---------------------------------------------------------------------- 
    4167   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    4672 
    4773   SUBROUTINE ice_thd_pnd 
     74 
    4875      !!------------------------------------------------------------------- 
    4976      !!               ***  ROUTINE ice_thd_pnd   *** 
    5077      !!                
    5178      !! ** Purpose :   change melt pond fraction and thickness 
    52       !!                 
     79      !! 
     80      !! ** Note    : Melt ponds affect only radiative transfer for now 
     81      !!              No heat, no salt. 
     82      !!              The current diagnostics lacks a contribution from drainage 
    5383      !!------------------------------------------------------------------- 
     84      INTEGER ::   ji, jj, jl        ! loop indices 
     85      !!------------------------------------------------------------------- 
     86       
     87      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 
     88      ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 
    5489      ! 
    55       SELECT CASE ( nice_pnd ) 
     90      diag_dvpn_mlt (:,:) = 0._wp   ;   diag_dvpn_drn (:,:) = 0._wp 
     91      diag_dvpn_lid (:,:) = 0._wp   ;   diag_dvpn_rnf (:,:) = 0._wp 
     92      diag_dvpn_mlt_1d(:) = 0._wp   ;   diag_dvpn_drn_1d(:) = 0._wp 
     93      diag_dvpn_lid_1d(:) = 0._wp   ;   diag_dvpn_rnf_1d(:) = 0._wp 
     94 
     95      !------------------------------------- 
     96      !  Remove ponds where ice has vanished 
     97      !------------------------------------- 
     98      at_i(:,:) = SUM( a_i, dim=3 ) 
    5699      ! 
    57       CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==! 
    58          ! 
    59       CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==! 
    60          ! 
    61       END SELECT 
     100      DO jl = 1, jpl 
     101         DO_2D( 1, 1, 1, 1 ) 
     102            IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN 
     103               wfx_pnd  (ji,jj)    = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 
     104               a_ip     (ji,jj,jl) = 0._wp 
     105               v_ip     (ji,jj,jl) = 0._wp 
     106               v_il     (ji,jj,jl) = 0._wp 
     107               h_ip     (ji,jj,jl) = 0._wp 
     108               h_il     (ji,jj,jl) = 0._wp 
     109               a_ip_frac(ji,jj,jl) = 0._wp 
     110            ENDIF 
     111         END_2D 
     112      END DO 
     113       
     114      !------------------------------ 
     115      !  Identify grid cells with ice 
     116      !------------------------------ 
     117      npti = 0   ;   nptidx(:) = 0 
     118      DO_2D( 1, 1, 1, 1 ) 
     119         IF( at_i(ji,jj) >= epsi10 ) THEN 
     120            npti = npti + 1 
     121            nptidx( npti ) = (jj - 1) * jpi + ji 
     122         ENDIF 
     123      END_2D 
     124 
     125      !------------------------------------ 
     126      !  Select melt pond scheme to be used 
     127      !------------------------------------ 
     128      IF( npti > 0 ) THEN 
     129         SELECT CASE ( nice_pnd ) 
     130            ! 
     131         CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==! 
     132            ! 
     133         CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
     134            ! 
     135         CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
     136            ! 
     137         END SELECT 
     138      ENDIF 
     139 
     140      !------------------------------------ 
     141      !  Diagnostics 
     142      !------------------------------------ 
     143      CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting 
     144      CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid 
     145      CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage 
     146      CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow 
    62147      ! 
     148      DEALLOCATE( diag_dvpn_mlt   , diag_dvpn_lid   , diag_dvpn_drn   , diag_dvpn_rnf    ) 
     149      DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 
     150       
    63151   END SUBROUTINE ice_thd_pnd  
    64152 
    … …  
    80168      !! ** References : Bush, G.W., and Trump, D.J. (2017) 
    81169      !!------------------------------------------------------------------- 
    82       INTEGER  ::   ji        ! loop indices 
     170      INTEGER  ::   ji, jl    ! loop indices 
     171      REAL(wp) ::   zdv_pnd   ! Amount of water going into the ponds & lids 
    83172      !!------------------------------------------------------------------- 
    84       DO ji = 1, npti 
    85          ! 
    86          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    87             h_ip_1d(ji)      = rn_hpnd     
    88             a_ip_1d(ji)      = rn_apnd * a_i_1d(ji) 
    89             h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever 
    90          ELSE 
    91             h_ip_1d(ji)      = 0._wp     
    92             a_ip_1d(ji)      = 0._wp 
    93             h_il_1d(ji)      = 0._wp 
    94          ENDIF 
    95          ! 
     173      DO jl = 1, jpl 
     174          
     175         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d    (1:npti), a_i    (:,:,jl) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d   (1:npti), t_su   (:,:,jl) ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) ) 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) ) 
     179         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) ) 
     180         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    ) 
     181 
     182         DO ji = 1, npti 
     183            ! 
     184            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
     185            ! 
     186            IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 
     187               h_ip_1d(ji)      = rn_hpnd     
     188               a_ip_1d(ji)      = rn_apnd * a_i_1d(ji) 
     189               h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever 
     190            ELSE 
     191               h_ip_1d(ji)      = 0._wp     
     192               a_ip_1d(ji)      = 0._wp 
     193               h_il_1d(ji)      = 0._wp 
     194            ENDIF 
     195            ! 
     196            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     197            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     198            ! 
     199            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     200            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 
     201            ! 
     202         END DO 
     203 
     204         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) ) 
     205         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) ) 
     206         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) ) 
     207         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d   (1:npti), v_ip   (:,:,jl) ) 
     208         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d   (1:npti), v_il   (:,:,jl) ) 
     209         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    ) 
     210 
    96211      END DO 
    97212      ! 
    … …  
    132247      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 --- 
    133248      !! 
    134       !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 --- 
    135       !!                                     perm = permability of sea-ice 
     249      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi * flush         --- from Flocco et al 2007 --- 
     250      !!                                     perm = permability of sea-ice                              + correction from Hunke et al 2012 (flush) 
    136251      !!                                     visc = water viscosity 
    137252      !!                                     Hp   = height of top of the pond above sea-level 
    138253      !!                                     Hi   = ice thickness thru which there is flushing 
     254      !!                                     flush= correction otherwise flushing is excessive 
    139255      !! 
    140256      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness 
    … …  
    143259      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect 
    144260      !! 
    145       !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 
     261      !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 
    146262      !!  
    147       !! ** Note       :   mostly stolen from CICE 
     263      !! ** Note       :   Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds.  
    148264      !! 
    149265      !! ** References :   Flocco and Feltham (JGR, 2007) 
    150266      !!                   Flocco et al       (JGR, 2010) 
    151267      !!                   Holland et al      (J. Clim, 2012) 
    152       !!------------------------------------------------------------------- 
     268      !!                   Hunke et al        (OM 2012) 
     269      !!-------------------------------------------------------------------   
    153270      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    154271      !! 
    … …  
    157274      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity 
    158275      !! 
    159       REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding 
     276      REAL(wp) ::   zfr_mlt, zdv_mlt, zdv_avail       ! fraction and volume of available meltwater retained for melt ponding 
    160277      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes 
     278      REAL(wp) ::   zdv_pnd                           ! Amount of water going into the ponds & lids 
    161279      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh 
    162280      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed 
    163281      REAL(wp) ::   zdT                               ! zTp-t_su 
    164       REAL(wp) ::   zsbr                              ! Brine salinity 
     282      REAL(wp) ::   zsbr, ztmelts                     ! Brine salinity 
    165283      REAL(wp) ::   zperm                             ! permeability of sea ice 
    166284      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    167285      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
    168286      !! 
    169       INTEGER  ::   ji, jk                            ! loop indices 
     287      INTEGER  ::   ji, jk, jl                        ! loop indices 
    170288      !!------------------------------------------------------------------- 
    171289      z1_rhow   = 1._wp / rhow  
    172290      z1_aspect = 1._wp / zaspect 
    173291      z1_Tp     = 1._wp / zTp  
    174  
    175       DO ji = 1, npti 
    176          !                                                            !----------------------------------------------------! 
    177          IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
    178             !                                                         !----------------------------------------------------! 
    179             !--- Remove ponds on thin ice or tiny ice fractions 
    180             a_ip_1d(ji)      = 0._wp 
    181             h_ip_1d(ji)      = 0._wp 
    182             h_il_1d(ji)      = 0._wp 
    183             !                                                         !--------------------------------! 
    184          ELSE                                                         ! Case ice thickness >= rn_himin ! 
    185             !                                                         !--------------------------------! 
    186             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     292       
     293      CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d          (1:npti), at_i          ) 
     294      CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d       (1:npti), wfx_pnd       ) 
     295       
     296      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     297      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     298      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     299      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
     300 
     301      DO jl = 1, jpl 
     302 
     303         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i (:,:,jl) ) 
     304         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d  (1:npti), h_i (:,:,jl) ) 
     305         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) ) 
     306         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) ) 
     307         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) ) 
     308         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) ) 
     309 
     310         CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) ) 
     311         CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) ) 
     312 
     313         DO jk = 1, nlay_i 
     314            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 
     315            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 
     316         END DO 
     317          
     318         !----------------------- 
     319         ! Melt pond calculations 
     320         !----------------------- 
     321         DO ji = 1, npti 
     322            ! 
     323            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
     324            !                                                            !----------------------------------------------------! 
     325            IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN   ! Case ice thickness < rn_himin or tiny ice fraction ! 
     326               !                                                         !----------------------------------------------------! 
     327               !--- Remove ponds on thin ice or tiny ice fractions 
     328               a_ip_1d(ji) = 0._wp 
     329               h_ip_1d(ji) = 0._wp 
     330               h_il_1d(ji) = 0._wp 
     331               !                                                         !--------------------------------! 
     332            ELSE                                                         ! Case ice thickness >= rn_himin ! 
     333               !                                                         !--------------------------------! 
     334               v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     335               v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     336               ! 
     337               !------------------! 
     338               ! case ice melting ! 
     339               !------------------! 
     340               ! 
     341               !--- available meltwater for melt ponding (zdv_avail) ---! 
     342               zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 
     343               zfr_mlt   = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
     344               zdv_mlt   = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors?  
     345               ! 
     346               !--- overflow ---! 
     347               ! 
     348               ! area driven overflow 
     349               !    If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     350               !       a_ip_max = zfr_mlt * a_i 
     351               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     352               zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     353               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     354 
     355               ! depth driven overflow 
     356               !    If pond depth exceeds half the ice thickness then reduce the pond volume 
     357               !       h_ip_max = 0.5 * h_i 
     358               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     359               zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
     360               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     361 
     362               !--- Pond growing ---! 
     363               v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     364               ! 
     365               !--- Lid melting ---! 
     366               IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     367               ! 
     368               !-------------------! 
     369               ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     370               !-------------------! 
     371               ! 
     372               zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
     373               ! 
     374               !--- Pond contraction (due to refreezing) ---! 
     375               IF( ln_pnd_lids ) THEN 
     376                  ! 
     377                  !--- Lid growing and subsequent pond shrinking ---!  
     378                  zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     379                     &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
     380 
     381                  ! Lid growing 
     382                  v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 
     383 
     384                  ! Pond shrinking 
     385                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
     386 
     387               ELSE 
     388                  zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp )  ! Holland 2012 (eq. 6)  
     389                  ! Pond shrinking 
     390                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
     391               ENDIF 
     392               ! 
     393               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     394               ! v_ip     = h_ip * a_ip 
     395               ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
     396               a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     397               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     398               ! 
     399 
     400               !------------------------------------------------!             
     401               ! Pond drainage through brine network (flushing) ! 
     402               !------------------------------------------------! 
     403               ! height of top of the pond above sea-level 
     404               zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
     405 
     406               ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
     407               DO jk = 1, nlay_i 
     408                  ! MV Assur is inconsistent with SI3 
     409                  !!zsbr = - 1.2_wp                                  & 
     410                  !!   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
     411                  !!   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
     412                  !!   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
     413                  !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     414                  ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
     415                  ztmelts  = -rTmlt * sz_i_1d(ji,jk) 
     416                  ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 
     417               END DO 
     418               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     419 
     420               ! Do the drainage using Darcy's law 
     421               zdv_flush   = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 
     422               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0  
     423               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     424 
     425               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     426               a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     427               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     428 
     429               !--- Corrections and lid thickness ---! 
     430               IF( ln_pnd_lids ) THEN 
     431                  !--- retrieve lid thickness from volume ---! 
     432                  IF( a_ip_1d(ji) > 0.01_wp ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     433                  ELSE                               ;   h_il_1d(ji) = 0._wp 
     434                  ENDIF 
     435                  !--- remove ponds if lids are much larger than ponds ---! 
     436                  IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     437                     a_ip_1d(ji) = 0._wp 
     438                     h_ip_1d(ji) = 0._wp 
     439                     h_il_1d(ji) = 0._wp 
     440                  ENDIF 
     441               ENDIF 
     442 
     443               ! diagnostics: dvpnd = mlt+rnf+lid+drn 
     444               diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow *   zdv_avail             * r1_Dt_ice   ! > 0, surface melt input 
     445               diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice   ! < 0, runoff 
     446               diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow *   zdv_frz               * r1_Dt_ice   ! < 0, shrinking 
     447               diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow *   zdv_flush             * r1_Dt_ice   ! < 0, drainage 
     448               ! 
     449            ENDIF 
     450            ! 
     451            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    187452            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
    188453            ! 
    189             !------------------! 
    190             ! case ice melting ! 
    191             !------------------! 
     454            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     455            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 
    192456            ! 
    193             !--- available meltwater for melt ponding ---! 
    194             zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    195             zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
    196             zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
    197             ! 
    198             !--- overflow ---! 
    199             ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    200             !    a_ip_max = zfr_mlt * a_i 
    201             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    202             zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    203             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    204  
    205             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    206             !    h_ip_max = 0.5 * h_i 
    207             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    208             zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    209             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    210              
    211             !--- Pond growing ---! 
    212             v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    213             ! 
    214             !--- Lid melting ---! 
    215             IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    216             ! 
    217             !--- mass flux ---! 
    218             IF( zdv_mlt > 0._wp ) THEN 
    219                zfac = zdv_mlt * rhow * r1_Dt_ice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    220                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    221                ! 
    222                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    223                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    224                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    225             ENDIF 
    226  
    227             !-------------------! 
    228             ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
    229             !-------------------! 
    230             ! 
    231             zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
    232             ! 
    233             !--- Pond contraction (due to refreezing) ---! 
    234             IF( ln_pnd_lids ) THEN 
    235                ! 
    236                !--- Lid growing and subsequent pond shrinking ---!  
    237                zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    238                   &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    239                 
    240                ! Lid growing 
    241                v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
    242                 
    243                ! Pond shrinking 
    244                v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    245  
    246             ELSE 
    247                ! Pond shrinking 
    248                v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
    249             ENDIF 
    250             ! 
    251             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    252             ! v_ip     = h_ip * a_ip 
    253             ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
    254             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    255             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    256  
    257             !---------------!             
    258             ! Pond flushing ! 
    259             !---------------! 
    260             ! height of top of the pond above sea-level 
    261             zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
    262              
    263             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
    264             DO jk = 1, nlay_i 
    265                zsbr = - 1.2_wp                                  & 
    266                   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
    267                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    268                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    269                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
    270             END DO 
    271             zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    272              
    273             ! Do the drainage using Darcy's law 
    274             zdv_flush   = -zperm * rho0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    275             zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    276             v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    277              
    278             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    279             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    280             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    281  
    282             !--- Corrections and lid thickness ---! 
    283             IF( ln_pnd_lids ) THEN 
    284                !--- retrieve lid thickness from volume ---! 
    285                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    286                ELSE                              ;   h_il_1d(ji) = 0._wp 
    287                ENDIF 
    288                !--- remove ponds if lids are much larger than ponds ---! 
    289                IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
    290                   a_ip_1d(ji)      = 0._wp 
    291                   h_ip_1d(ji)      = 0._wp 
    292                   h_il_1d(ji)      = 0._wp 
    293                ENDIF 
    294             ENDIF 
    295             ! 
    296          ENDIF 
    297           
     457         END DO 
     458 
     459         !-------------------------------------------------------------------- 
     460         ! Retrieve 2D arrays 
     461         !-------------------------------------------------------------------- 
     462         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 
     463         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 
     464         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 
     465         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 
     466         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 
     467         ! 
    298468      END DO 
    299469      ! 
     470      CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 
     471      ! 
     472      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     473      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     474      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     475      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        ) 
     476      ! 
    300477   END SUBROUTINE pnd_LEV 
    301478 
     479 
     480 
     481   SUBROUTINE pnd_TOPO     
     482                                          
     483      !!------------------------------------------------------------------- 
     484      !!                ***  ROUTINE pnd_TOPO  *** 
     485      !! 
     486      !! ** Purpose :   Compute melt pond evolution based on the ice 
     487      !!                topography inferred from the ice thickness distribution 
     488      !! 
     489      !! ** Method  :   This code is initially based on Flocco and Feltham 
     490      !!                (2007) and Flocco et al. (2010).  
     491      !! 
     492      !!                - Calculate available pond water base on surface meltwater 
     493      !!                - Redistribute water as a function of topography, drain water 
     494      !!                - Exchange water with the lid 
     495      !! 
     496      !! ** Tunable parameters : 
     497      !! 
     498      !! ** Note : 
     499      !! 
     500      !! ** References 
     501      !! 
     502      !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond 
     503      !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi: 
     504      !!    10.1029/2006JC003836. 
     505      !! 
     506      !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
     507      !!    a physically based melt pond scheme into the sea ice component of a 
     508      !!    climate model.  J. Geophys. Res. 115, C08012, 
     509      !!    doi: 10.1029/2009JC005568. 
     510      !! 
     511      !!------------------------------------------------------------------- 
     512      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     513         zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C) 
     514         zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     515 
     516 
     517      ! local variables 
     518      REAL(wp) :: & 
     519         zdHui,   &      ! change in thickness of ice lid (m) 
     520         zomega,  &      ! conduction 
     521         zdTice,  &      ! temperature difference across ice lid (C) 
     522         zdvice,  &      ! change in ice volume (m) 
     523         zTavg,   &      ! mean surface temperature across categories (C) 
     524         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2) 
     525         zTp,     &      ! pond freezing temperature (C) 
     526         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3) 
     527         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding 
     528         z1_rhow, &      ! inverse water density 
     529         zv_pnd  , &     ! volume of meltwater contributing to ponds 
     530         zv_mlt          ! total amount of meltwater produced 
     531 
     532      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, &     !! total melt pond water available before redistribution and drainage 
     533                                        zvolp_res    !! remaining melt pond water available after drainage 
     534                                         
     535      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i 
     536 
     537      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     538 
     539      INTEGER  ::   i_test 
     540 
     541      ! Note 
     542      ! equivalent for CICE translation 
     543      ! a_ip      -> apond 
     544      ! a_ip_frac -> apnd 
     545 
     546      CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 
     547       
     548      !--------------------------------------------------------------- 
     549      ! Initialise 
     550      !--------------------------------------------------------------- 
     551 
     552      ! Parameters & constants (move to parameters) 
     553      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3) 
     554      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline) 
     555      z1_rhow   = 1._wp / rhow  
     556 
     557      ! Set required ice variables (hard-coded here for now) 
     558!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?)  
     559       
     560      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 
     561      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 
     562      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 
     563      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 
     564       
     565      ! thickness 
     566      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 
     567      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp 
     568      END WHERE 
     569      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 
     570       
     571      !--------------------------------------------------------------- 
     572      ! Change 2D to 1D 
     573      !--------------------------------------------------------------- 
     574      ! MV  
     575      ! a less computing-intensive version would have 2D-1D passage here 
     576      ! use what we have in iceitd.F90 (incremental remapping) 
     577 
     578      !-------------------------------------------------------------- 
     579      ! Collect total available pond water volume 
     580      !-------------------------------------------------------------- 
     581      ! Assuming that meltwater (+rain in principle) runsoff the surface 
     582      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 
     583      ! I cite her words, they are very talkative 
     584      ! "grid cells with very little ice cover (and hence more open water area)  
     585      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 
     586      ! "This results in the same runoff fraction r for each ice category within a grid cell" 
     587       
     588      zvolp(:,:) = 0._wp 
     589 
     590      DO jl = 1, jpl 
     591         DO_2D( 1, 1, 1, 1 ) 
     592                  
     593               IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
     594             
     595                  !--- Available and contributing meltwater for melt ponding ---! 
     596                  zv_mlt  = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) &        ! available volume of surface melt water per grid area 
     597                     &    * z1_rhow * a_i(ji,jj,jl) 
     598                      ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 
     599                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds 
     600                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl 
     601 
     602                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice                     ! diags 
     603                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice    
     604 
     605                  !--- Create possible new ponds 
     606                  ! if pond does not exist, create new pond over full ice area 
     607                  !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 
     608                  IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 
     609                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl) 
     610                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE) 
     611                  ENDIF 
     612                   
     613                  !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 
     614                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness 
     615                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     616                   
     617                  !--- Total available pond water volume (pre-existing + newly produced)j 
     618                  zvolp(ji,jj)   = zvolp(ji,jj)   + v_ip(ji,jj,jl)  
     619!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 
     620                    
     621               ENDIF ! a_i 
     622 
     623         END_2D 
     624      END DO ! ji 
     625                   
     626      !-------------------------------------------------------------- 
     627      ! Redistribute and drain water from ponds 
     628      !--------------------------------------------------------------    
     629      CALL ice_thd_pnd_area( zvolp, zvolp_res ) 
     630                                    
     631      !-------------------------------------------------------------- 
     632      ! Melt pond lid growth and melt 
     633      !--------------------------------------------------------------    
     634       
     635      IF( ln_pnd_lids ) THEN 
     636 
     637         DO_2D( 1, 1, 1, 1 ) 
     638 
     639            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     640                   
     641               !-------------------------- 
     642               ! Pond lid growth and melt 
     643               !-------------------------- 
     644               ! Mean surface temperature 
     645               zTavg = 0._wp 
     646               DO jl = 1, jpl 
     647                  zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 
     648               END DO 
     649               zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 
     650          
     651               DO jl = 1, jpl-1 
     652             
     653                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN 
     654                
     655                     !---------------------------------------------------------------- 
     656                     ! Lid melting: floating upper ice layer melts in whole or part 
     657                     !---------------------------------------------------------------- 
     658                     ! Use Tsfc for each category 
     659                     IF ( t_su(ji,jj,jl) > zTp ) THEN 
     660 
     661                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 
     662                         
     663                        IF ( zdvice > epsi10 ) THEN 
     664                         
     665                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice 
     666                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -  
     667                                                                        ! as it is already counted in surface melt 
     668!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar) 
     669!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw) 
     670                      
     671                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 
     672                           ! ice lid melted and category is pond covered 
     673                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl)  
     674!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl)  
     675                              v_il(ji,jj,jl)   = 0._wp 
     676                           ENDIF 
     677                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 
     678                            
     679                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag 
     680                            
     681                        ENDIF 
     682                         
     683                     !---------------------------------------------------------------- 
     684                     ! Freeze pre-existing lid  
     685                     !---------------------------------------------------------------- 
     686 
     687                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 
     688 
     689                        ! differential growth of base of surface floating ice layer 
     690                        zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0    
     691                        zomega = rcnd_i * zdTice / zrhoi_L 
     692                        zdHui  = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 
     693                               - v_il(ji,jj,jl) / a_i(ji,jj,jl) 
     694                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     695                   
     696                        IF ( zdvice > epsi10 ) THEN 
     697                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     698                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     699!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice 
     700!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice 
     701                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     702                            
     703                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag 
     704                            
     705                        ENDIF 
     706                   
     707                     ENDIF ! Tsfcn(i,j,n) 
     708 
     709                  !---------------------------------------------------------------- 
     710                  ! Freeze new lids 
     711                  !---------------------------------------------------------------- 
     712                  !  upper ice layer begins to form 
     713                  ! note: albedo does not change 
     714 
     715                  ELSE ! v_il < epsi10 
     716                     
     717                     ! thickness of newly formed ice 
     718                     ! the surface temperature of a meltpond is the same as that 
     719                     ! of the ice underneath (0C), and the thermodynamic surface  
     720                     ! flux is the same 
     721                      
     722                     !!! we need net surface energy flux, excluding conduction 
     723                     !!! fsurf is summed over categories in CICE 
     724                     !!! we have the category-dependent flux, let us use it ? 
     725                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                      
     726                     zdHui  = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) 
     727                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     728                     IF ( zdvice > epsi10 ) THEN 
     729                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     730                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     731                         
     732                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag 
     733!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice 
     734!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice 
     735                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 
     736                     ENDIF 
     737                
     738                  ENDIF  ! v_il 
     739             
     740               END DO ! jl 
     741 
     742            ELSE  ! remove ponds on thin ice 
     743 
     744               v_ip(ji,jj,:) = 0._wp 
     745               v_il(ji,jj,:) = 0._wp 
     746!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 
     747!                 zvolp(ji,jj)    = 0._wp          
     748 
     749            ENDIF 
     750 
     751         END_2D 
     752 
     753      ENDIF ! ln_pnd_lids 
     754 
     755      !--------------------------------------------------------------- 
     756      ! Clean-up variables (probably duplicates what icevar would do) 
     757      !--------------------------------------------------------------- 
     758      ! MV comment 
     759      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 
     760      ! icevar should recompute all other variables (if needed at all) 
     761 
     762      DO jl = 1, jpl 
     763 
     764         DO_2D( 1, 1, 1, 1 ) 
     765 
     766!              ! zap lids on small ponds 
     767!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 
     768!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN 
     769!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 
     770!              ENDIF 
     771       
     772               ! recalculate equivalent pond variables 
     773               IF ( a_ip(ji,jj,jl) > epsi10) THEN 
     774                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 
     775                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     776                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     777               ENDIF 
     778!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar 
     779!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar 
     780!              ENDIF 
     781                
     782         END_2D 
     783 
     784      END DO   ! jl 
     785 
     786 
     787   END SUBROUTINE pnd_TOPO 
     788 
     789    
     790   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 
     791 
     792       !!------------------------------------------------------------------- 
     793       !!                ***  ROUTINE ice_thd_pnd_area *** 
     794       !! 
     795       !! ** Purpose : Given the total volume of available pond water,  
     796       !!              redistribute and drain water 
     797       !! 
     798       !! ** Method 
     799       !! 
     800       !-----------| 
     801       !           | 
     802       !           |-----------| 
     803       !___________|___________|______________________________________sea-level 
     804       !           |           | 
     805       !           |           |---^--------| 
     806       !           |           |   |        | 
     807       !           |           |   |        |-----------|              |------- 
     808       !           |           |   | alfan  |           |              | 
     809       !           |           |   |        |           |--------------| 
     810       !           |           |   |        |           |              | 
     811       !---------------------------v------------------------------------------- 
     812       !           |           |   ^        |           |              | 
     813       !           |           |   |        |           |--------------| 
     814       !           |           |   | betan  |           |              | 
     815       !           |           |   |        |-----------|              |------- 
     816       !           |           |   |        | 
     817       !           |           |---v------- | 
     818       !           |           | 
     819       !           |-----------| 
     820       !           | 
     821       !-----------| 
     822       ! 
     823       !! 
     824       !!------------------------------------------------------------------ 
     825        
     826       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 
     827          zvolp,                                       &  ! total available pond water 
     828          zdvolp                                          ! remaining meltwater after redistribution 
     829 
     830       INTEGER ::  & 
     831          ns,      & 
     832          m_index, & 
     833          permflag 
     834 
     835       REAL (wp), DIMENSION(jpl) :: & 
     836          hicen, & 
     837          hsnon, & 
     838          asnon, & 
     839          alfan, & 
     840          betan, & 
     841          cum_max_vol, & 
     842          reduced_aicen 
     843 
     844       REAL (wp), DIMENSION(0:jpl) :: & 
     845          cum_max_vol_tmp 
     846 
     847       REAL (wp) :: & 
     848          hpond, & 
     849          drain, & 
     850          floe_weight, & 
     851          pressure_head, & 
     852          hsl_rel, & 
     853          deltah, & 
     854          perm, & 
     855          msno 
     856 
     857       REAL (wp), parameter :: & 
     858          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s 
     859 
     860      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     861         zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     862 
     863      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     864 
     865       a_ip(:,:,:) = 0._wp 
     866       h_ip(:,:,:) = 0._wp 
     867  
     868       DO_2D( 1, 1, 1, 1 ) 
     869  
     870             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     871  
     872        !------------------------------------------------------------------- 
     873        ! initialize 
     874        !------------------------------------------------------------------- 
     875  
     876        DO jl = 1, jpl 
     877  
     878           !---------------------------------------- 
     879           ! compute the effective snow fraction 
     880           !---------------------------------------- 
     881  
     882           IF (a_i(ji,jj,jl) < epsi10)  THEN 
     883              hicen(jl) =  0._wp 
     884              hsnon(jl) =  0._wp 
     885              reduced_aicen(jl) = 0._wp 
     886              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 
     887           ELSE 
     888              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     889              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     890              reduced_aicen(jl) = 1._wp ! n=jpl 
     891  
     892              !js: initial code in NEMO_DEV 
     893              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 
     894              !                     * (-0.024_wp*hicen(jl) + 0.832_wp) 
     895  
     896              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 
     897              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &  
     898                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 
     899  
     900              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical) 
     901              ! MV should check whether this makes sense to have the same effective snow fraction in here 
     902              ! OLI: it probably doesn't 
     903           END IF 
     904  
     905  ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     906  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     907  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     908  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     909  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     910  ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
     911  ! alfan = 60% of the ice volume) in each category lies above the reference line, 
     912  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     913  
     914  ! MV: 
     915  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
     916  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
     917  
     918  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 
     919  
     920           alfan(jl) = 0.6 * hicen(jl) 
     921           betan(jl) = 0.4 * hicen(jl) 
     922  
     923           cum_max_vol(jl)     = 0._wp 
     924           cum_max_vol_tmp(jl) = 0._wp 
     925  
     926        END DO ! jpl 
     927  
     928        cum_max_vol_tmp(0) = 0._wp 
     929        drain = 0._wp 
     930        zdvolp(ji,jj) = 0._wp 
     931  
     932        !---------------------------------------------------------- 
     933        ! Drain overflow water, update pond fraction and volume 
     934        !---------------------------------------------------------- 
     935  
     936        !-------------------------------------------------------------------------- 
     937        ! the maximum amount of water that can be contained up to each ice category 
     938        !-------------------------------------------------------------------------- 
     939        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     940        ! Then the excess volume cum_max_vol(jl) drains out of the system 
     941        ! It should be added to wfx_pnd_out 
     942  
     943        DO jl = 1, jpl-1 ! last category can not hold any volume 
     944  
     945           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 
     946  
     947              ! total volume in level including snow 
     948              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 
     949                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 
     950  
     951              ! subtract snow solid volumes from lower categories in current level 
     952              DO ns = 1, jl 
     953                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 
     954                    - rhos/rhow * &     ! free air fraction that can be filled by water 
     955                      asnon(ns)  * &    ! effective areal fraction of snow in that category 
     956                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 
     957              END DO 
     958  
     959           ELSE ! assume higher categories unoccupied 
     960              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 
     961           END IF 
     962           !IF (cum_max_vol_tmp(jl) < z0) THEN 
     963           !   CALL abort_ice('negative melt pond volume') 
     964           !END IF 
     965        END DO 
     966        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     967        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     968  
     969        !---------------------------------------------------------------- 
     970        ! is there more meltwater than can be held in the floe? 
     971        !---------------------------------------------------------------- 
     972        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 
     973           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 
     974           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 
     975  
     976           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice) 
     977            
     978           zdvolp(ji,jj) = drain         ! this is the drained water 
     979           IF (zvolp(ji,jj) < epsi10) THEN 
     980              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     981              zvolp(ji,jj) = 0._wp 
     982           END IF 
     983        END IF 
     984  
     985        ! height and area corresponding to the remaining volume 
     986        ! routine leaves zvolp unchanged 
     987        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     988  
     989        DO jl = 1, m_index 
     990           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 
     991           !                                         !  volume instead, no ? 
     992           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2 
     993           a_ip(ji,jj,jl) = reduced_aicen(jl) 
     994           ! in practise, pond fraction depends on the empirical snow fraction 
     995           ! so in turn on ice thickness 
     996        END DO 
     997        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     998  
     999        !------------------------------------------------------------------------ 
     1000        ! Drainage through brine network (permeability) 
     1001        !------------------------------------------------------------------------ 
     1002        !!! drainage due to ice permeability - Darcy's law 
     1003  
     1004        ! sea water level 
     1005        msno = 0._wp  
     1006        DO jl = 1 , jpl 
     1007          msno = msno + v_s(ji,jj,jl) * rhos 
     1008        END DO 
     1009        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj) 
     1010        hsl_rel = floe_weight / rho0 & 
     1011                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 
     1012  
     1013        deltah = hpond - hsl_rel 
     1014        pressure_head = grav * rho0 * max(deltah, 0._wp) 
     1015  
     1016        ! drain if ice is permeable 
     1017        permflag = 0 
     1018  
     1019        IF (pressure_head > 0._wp) THEN 
     1020           DO jl = 1, jpl-1 
     1021              IF ( hicen(jl) /= 0._wp ) THEN 
     1022  
     1023              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2 
     1024  
     1025                 perm = 0._wp ! MV ugly dummy patch 
     1026                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
     1027                 IF (perm > 0._wp) permflag = 1 
     1028  
     1029                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & 
     1030                                          (viscosity*hicen(jl)) 
     1031                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 
     1032                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 
     1033  
     1034                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 
     1035                  
     1036                 IF (zvolp(ji,jj) < epsi10) THEN 
     1037                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     1038                    zvolp(ji,jj) = 0._wp 
     1039                 END IF 
     1040             END IF 
     1041          END DO 
     1042  
     1043           ! adjust melt pond dimensions 
     1044           IF (permflag > 0) THEN 
     1045              ! recompute pond depth 
     1046             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     1047              DO jl = 1, m_index 
     1048                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 
     1049                 a_ip(ji,jj,jl) = reduced_aicen(jl) 
     1050              END DO 
     1051              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1052           END IF 
     1053        END IF ! pressure_head 
     1054  
     1055        !------------------------------- 
     1056        ! remove water from the snow 
     1057        !------------------------------- 
     1058        !------------------------------------------------------------------------ 
     1059        ! total melt pond volume in category does not include snow volume 
     1060        ! snow in melt ponds is not melted 
     1061        !------------------------------------------------------------------------ 
     1062         
     1063        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 
     1064        ! how much, so I did not diagnose it 
     1065        ! so if there is a problem here, nobody is going to see it... 
     1066         
     1067  
     1068        ! Calculate pond volume for lower categories 
     1069        DO jl = 1,m_index-1 
     1070           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 
     1071                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
     1072        END DO 
     1073  
     1074        ! Calculate pond volume for highest category = remaining pond volume 
     1075  
     1076        ! The following is completely unclear to Martin at least 
     1077        ! Could we redefine properly and recode in a more readable way ? 
     1078  
     1079        ! m_index = last category with melt pond 
     1080  
     1081        IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 
     1082  
     1083        IF (m_index > 1) THEN 
     1084          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 
     1085             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 
     1086          ELSE 
     1087             v_ip(ji,jj,m_index) = 0._wp  
     1088             h_ip(ji,jj,m_index) = 0._wp 
     1089             a_ip(ji,jj,m_index) = 0._wp 
     1090             ! If remaining pond volume is negative reduce pond volume of 
     1091             ! lower category 
     1092             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 
     1093              v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 
     1094          END IF 
     1095        END IF 
     1096  
     1097        DO jl = 1,m_index 
     1098           IF (a_ip(ji,jj,jl) > epsi10) THEN 
     1099               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     1100           ELSE 
     1101              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 
     1102              h_ip(ji,jj,jl) = 0._wp  
     1103              v_ip(ji,jj,jl)  = 0._wp 
     1104              a_ip(ji,jj,jl) = 0._wp 
     1105           END IF 
     1106        END DO 
     1107        DO jl = m_index+1, jpl 
     1108           h_ip(ji,jj,jl) = 0._wp  
     1109           a_ip(ji,jj,jl) = 0._wp  
     1110           v_ip(ji,jj,jl) = 0._wp  
     1111        END DO 
     1112         
     1113           ENDIF 
     1114 
     1115     END_2D 
     1116 
     1117    END SUBROUTINE ice_thd_pnd_area 
     1118 
     1119 
     1120    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 
     1121       !!------------------------------------------------------------------- 
     1122       !!                ***  ROUTINE ice_thd_pnd_depth  *** 
     1123       !! 
     1124       !! ** Purpose :   Compute melt pond depth 
     1125       !!------------------------------------------------------------------- 
     1126 
     1127       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     1128          aicen, & 
     1129          asnon, & 
     1130          hsnon, & 
     1131          alfan, & 
     1132          cum_max_vol 
     1133 
     1134       REAL (wp), INTENT(IN) :: & 
     1135          zvolp 
     1136 
     1137       REAL (wp), INTENT(OUT) :: & 
     1138          hpond 
     1139 
     1140       INTEGER, INTENT(OUT) :: & 
     1141          m_index 
     1142 
     1143       INTEGER :: n, ns 
     1144 
     1145       REAL (wp), DIMENSION(0:jpl+1) :: & 
     1146          hitl, & 
     1147          aicetl 
     1148 
     1149       REAL (wp) :: & 
     1150          rem_vol, & 
     1151          area, & 
     1152          vol, & 
     1153          tmp, & 
     1154          z0   = 0.0_wp 
     1155 
     1156       !---------------------------------------------------------------- 
     1157       ! hpond is zero if zvolp is zero - have we fully drained? 
     1158       !---------------------------------------------------------------- 
     1159 
     1160       IF (zvolp < epsi10) THEN 
     1161        hpond = z0 
     1162        m_index = 0 
     1163       ELSE 
     1164 
     1165        !---------------------------------------------------------------- 
     1166        ! Calculate the category where water fills up to 
     1167        !---------------------------------------------------------------- 
     1168 
     1169        !----------| 
     1170        !          | 
     1171        !          | 
     1172        !          |----------|                                     -- -- 
     1173        !__________|__________|_________________________________________ ^ 
     1174        !          |          |             rem_vol     ^                | Semi-filled 
     1175        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     1176        !          |          |          |              | 
     1177        !          |          |          |              |hpond 
     1178        !          |          |          |----------|   |     |------- 
     1179        !          |          |          |          |   |     | 
     1180        !          |          |          |          |---v-----| 
     1181        !          |          | m_index  |          |         | 
     1182        !------------------------------------------------------------- 
     1183 
     1184        m_index = 0  ! 1:m_index categories have water in them 
     1185        DO n = 1, jpl 
     1186           IF (zvolp <= cum_max_vol(n)) THEN 
     1187              m_index = n 
     1188              IF (n == 1) THEN 
     1189                 rem_vol = zvolp 
     1190              ELSE 
     1191                 rem_vol = zvolp - cum_max_vol(n-1) 
     1192              END IF 
     1193              exit ! to break out of the loop 
     1194           END IF 
     1195        END DO 
     1196        m_index = min(jpl-1, m_index) 
     1197 
     1198        !---------------------------------------------------------------- 
     1199        ! semi-filled layer may have m_index different snow in it 
     1200        !---------------------------------------------------------------- 
     1201 
     1202        !-----------------------------------------------------------  ^ 
     1203        !                                                             |  alfan(m_index+1) 
     1204        !                                                             | 
     1205        !hitl(3)-->                             |----------|          | 
     1206        !hitl(2)-->                |------------| * * * * *|          | 
     1207        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     1208        !hitl(0)-->-------------------------------------------------  |  ^ 
     1209        !                various snow from lower categories          |  |alfa(m_index) 
     1210 
     1211        ! hitl - heights of the snow layers from thinner and current categories 
     1212        ! aicetl - area of each snow depth in this layer 
     1213 
     1214        hitl(:) = z0 
     1215        aicetl(:) = z0 
     1216        DO n = 1, m_index 
     1217           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     1218                                  alfan(m_index+1) - alfan(m_index)), z0) 
     1219           aicetl(n) = asnon(n) 
     1220 
     1221           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     1222        END DO 
     1223 
     1224        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     1225        aicetl(m_index+1) = z0 
     1226 
     1227        !---------------------------------------------------------------- 
     1228        ! reorder array according to hitl 
     1229        ! snow heights not necessarily in height order 
     1230        !---------------------------------------------------------------- 
     1231 
     1232        DO ns = 1, m_index+1 
     1233           DO n = 0, m_index - ns + 1 
     1234              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     1235                 tmp = hitl(n) 
     1236                 hitl(n) = hitl(n+1) 
     1237                 hitl(n+1) = tmp 
     1238                 tmp = aicetl(n) 
     1239                 aicetl(n) = aicetl(n+1) 
     1240                 aicetl(n+1) = tmp 
     1241              END IF 
     1242           END DO 
     1243        END DO 
     1244 
     1245        !---------------------------------------------------------------- 
     1246        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     1247        !---------------------------------------------------------------- 
     1248 
     1249        !hitl(3)---------------------------------------------------------------- 
     1250        !                                                       | * * * * * * * * 
     1251        !                                                       |* * * * * * * * * 
     1252        !hitl(2)---------------------------------------------------------------- 
     1253        !                                    | * * * * * * * *  | * * * * * * * * 
     1254        !                                    |* * * * * * * * * |* * * * * * * * * 
     1255        !hitl(1)---------------------------------------------------------------- 
     1256        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
     1257        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
     1258        !hitl(0)---------------------------------------------------------------- 
     1259        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3) 
     1260 
     1261        ! move up over layers incrementing volume 
     1262        DO n = 1, m_index+1 
     1263 
     1264           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     1265                (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     1266 
     1267           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     1268 
     1269           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     1270              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 
     1271 
     1272              exit 
     1273           ELSE  ! still in sub-layer below the sub-layer with the depth 
     1274              rem_vol = rem_vol - vol 
     1275           END IF 
     1276 
     1277        END DO 
     1278 
     1279       END IF 
     1280 
     1281    END SUBROUTINE ice_thd_pnd_depth 
     1282 
     1283 
     1284    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 
     1285       !!------------------------------------------------------------------- 
     1286       !!                ***  ROUTINE ice_thd_pnd_perm *** 
     1287       !! 
     1288       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     1289       !!                and its permeability 
     1290       !!------------------------------------------------------------------- 
     1291 
     1292       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     1293          ticen,  &     ! internal ice temperature (K) 
     1294          salin         ! salinity (ppt)     !js: ppt according to cice 
     1295 
     1296       REAL (wp), INTENT(OUT) :: & 
     1297          perm      ! permeability 
     1298 
     1299       REAL (wp) ::   & 
     1300          Sbr       ! brine salinity 
     1301 
     1302       REAL (wp), DIMENSION(nlay_i) ::   & 
     1303          Tin, &    ! ice temperature 
     1304          phi       ! liquid fraction 
     1305 
     1306       INTEGER :: k 
     1307 
     1308       !----------------------------------------------------------------- 
     1309       ! Compute ice temperatures from enthalpies using quadratic formula 
     1310       !----------------------------------------------------------------- 
     1311 
     1312       DO k = 1,nlay_i 
     1313          Tin(k) = ticen(k) - rt0   !js: from K to degC 
     1314       END DO 
     1315 
     1316       !----------------------------------------------------------------- 
     1317       ! brine salinity and liquid fraction 
     1318       !----------------------------------------------------------------- 
     1319 
     1320       DO k = 1, nlay_i 
     1321        
     1322          Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 
     1323          ! Best expression to date is that one (Vancoppenolle et al JGR 2019) 
     1324          ! Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 
     1325          phi(k) = salin(k) / Sbr 
     1326           
     1327       END DO 
     1328 
     1329       !----------------------------------------------------------------- 
     1330       ! permeability 
     1331       !----------------------------------------------------------------- 
     1332 
     1333       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 
     1334 
     1335   END SUBROUTINE ice_thd_pnd_perm 
    3021336 
    3031337   SUBROUTINE ice_thd_pnd_init  
    … …  
    3151349      INTEGER  ::   ios, ioptio   ! Local integer 
    3161350      !! 
    317       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 
     1351      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 
    3181352         &                          ln_pnd_CST , rn_apnd, rn_hpnd,         & 
     1353         &                          ln_pnd_TOPO,                           & 
    3191354         &                          ln_pnd_lids, ln_pnd_alb 
    3201355      !!------------------------------------------------------------------- 
    … …  
    3321367         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    3331368         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd 
     1369         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO 
    3341370         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV 
    3351371         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min 
    3361372         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max 
     1373         WRITE(numout,*) '            Pond flushing efficiency                              rn_pnd_flush = ', rn_pnd_flush 
    3371374         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST 
    3381375         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd 
    … …  
    3471384      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    3481385      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF 
     1386      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF 
    3491387      IF( ioptio /= 1 )   & 
    350          & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' ) 
     1388         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' ) 
    3511389      ! 
    3521390      SELECT CASE( nice_pnd ) 

  • NEMO/trunk/src/ICE/icethd_zdf_bl99.F90

    r13472 r14005  
    109109      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function 
    110110      ! 
    111       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    112       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    113       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    114       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
    115       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    116       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    117       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    119       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    120       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    121       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    122       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    123       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    124       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    125       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    126       REAL(wp), DIMENSION(jpij)            ::   zkappa_comb ! Combined snow and ice surface conductivity 
    127       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    128       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    129       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    130       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    131       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    132       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    133       REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow  
    134       REAL(wp), DIMENSION(jpij)            ::   isnow       ! snow presence (1) or not (0)  
    135       REAL(wp), DIMENSION(jpij)            ::   isnow_comb  ! snow presence for met-office  
     111      REAL(wp), DIMENSION(jpij       )   ::   ztsuold     ! Old surface temperature in the ice 
     112      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztiold      ! Old temperature in the ice 
     113      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsold      ! Old temperature in the snow 
     114      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Temporary temperature in the ice to check the convergence 
     115      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow to check the convergence 
     116      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
     117      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i_cp ! copy 
     118      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
     119      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
     120      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
     121      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
     122      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
     123      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
     124      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
     125      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zeta_s      ! Eta factor in the snow 
     126      REAL(wp), DIMENSION(jpij)          ::   zkappa_comb ! Combined snow and ice surface conductivity 
     127      REAL(wp), DIMENSION(jpij)          ::   zq_ini      ! diag errors on heat 
     128      REAL(wp), DIMENSION(jpij)          ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     129      REAL(wp), DIMENSION(jpij)          ::   za_s_fra    ! ice fraction covered by snow  
     130      REAL(wp), DIMENSION(jpij)          ::   isnow       ! snow presence (1) or not (0)  
     131      REAL(wp), DIMENSION(jpij)          ::   isnow_comb  ! snow presence for met-office  
     132      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindterm    ! 'Ind'ependent term 
     133      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindtbis    ! Temporary 'ind'ependent term 
     134      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     135      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) ::   ztrid       ! Tridiagonal system terms 
    136136      ! 
    137137      ! Mono-category 
    … …  
    533533            ! Solve the tridiagonal system with Gauss elimination method. 
    534534            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    535             jm_maxt = 0 
    536             jm_mint = nlay_i+5 
    537             DO ji = 1, npti 
    538                jm_mint = MIN(jm_min(ji),jm_mint) 
    539                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    540             END DO 
    541  
    542             DO jk = jm_mint+1, jm_maxt 
    543                DO ji = 1, npti 
    544                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     535!!$            jm_maxt = 0 
     536!!$            jm_mint = nlay_i+5 
     537!!$            DO ji = 1, npti 
     538!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     539!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     540!!$            END DO 
     541!!$            !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 
     542!!$             
     543!!$            DO jk = jm_mint+1, jm_maxt 
     544!!$               DO ji = 1, npti 
     545!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     546!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     547!!$                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     548!!$               END DO 
     549!!$            END DO 
     550            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     551            DO ji = 1, npti 
     552               jm_mint = jm_min(ji) 
     553               jm_maxt = jm_max(ji) 
     554               DO jm = jm_mint+1, jm_maxt 
    545555                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    546556                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
    … …  
    564574            END DO 
    565575 
     576            ! snow temperatures       
    566577            DO ji = 1, npti 
    567578               ! Variables used after iterations 
    568579               ! Value must be frozen after convergence for MPP independance reason 
    569                IF ( .NOT. l_T_converged(ji) ) THEN 
    570                   ! snow temperatures       
    571                   IF( h_s_1d(ji) > 0._wp ) THEN 
    572                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    573                   ENDIF 
    574                   ! surface temperature 
     580               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     581                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     582            END DO 
     583            !!clem SNWLAY 
     584            DO jm = nlay_s, 2, -1 
     585               DO ji = 1, npti 
     586                  jk = jm - 1 
     587                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     588                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     589               END DO 
     590            END DO 
     591             
     592            ! surface temperature 
     593            DO ji = 1, npti 
     594               IF( .NOT. l_T_converged(ji) ) THEN 
    575595                  ztsub(ji) = t_su_1d(ji) 
    576596                  IF( t_su_1d(ji) < rt0 ) THEN 
    577                      t_su_1d(ji) = (  zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    578                         &           ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     597                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     598                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    579599                  ENDIF 
    580600               ENDIF 
    581601            END DO 
    582             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
    583602            ! 
    584603            !-------------------------------------------------------------- 
    … …  
    727746            ! Solve the tridiagonal system with Gauss elimination method. 
    728747            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    729             jm_maxt = 0 
    730             jm_mint = nlay_i+5 
    731             DO ji = 1, npti 
    732                jm_mint = MIN(jm_min(ji),jm_mint) 
    733                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    734             END DO 
    735              
    736             DO jk = jm_mint+1, jm_maxt 
    737                DO ji = 1, npti 
    738                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     748!!$            jm_maxt = 0 
     749!!$            jm_mint = nlay_i+5 
     750!!$            DO ji = 1, npti 
     751!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     752!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     753!!$            END DO 
     754!!$             
     755!!$            DO jk = jm_mint+1, jm_maxt 
     756!!$               DO ji = 1, npti 
     757!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     758!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     759!!$                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1) 
     760!!$               END DO 
     761!!$            END DO 
     762            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     763            DO ji = 1, npti 
     764               jm_mint = jm_min(ji) 
     765               jm_maxt = jm_max(ji) 
     766               DO jm = jm_mint+1, jm_maxt 
    739767                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    740                   zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)  / zdiagbis(ji,jm-1) 
     768                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
    741769               END DO 
    742770            END DO 
    743              
     771 
    744772            ! ice temperatures 
    745773            DO ji = 1, npti 
    … …  
    761789            ! snow temperatures       
    762790            DO ji = 1, npti 
    763                ! Variable used after iterations 
     791               ! Variables used after iterations 
    764792               ! Value must be frozen after convergence for MPP independance reason 
    765                IF ( .NOT. l_T_converged(ji) ) THEN 
    766                   IF( h_s_1d(ji) > 0._wp ) THEN 
    767                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    768                   ENDIF 
    769                ENDIF 
    770             END DO 
    771             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
     793               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     794                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     795            END DO 
     796            !!clem SNWLAY 
     797            DO jm = nlay_s, 2, -1 
     798               DO ji = 1, npti 
     799                  jk = jm - 1 
     800                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     801                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     802               END DO 
     803            END DO 
    772804            ! 
    773805            !-------------------------------------------------------------- 
    … …  
    923955         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    924956         IF( h_s_1d(ji) >= zhs_ssl ) THEN 
    925             t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1)   & 
    926                &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 
     957            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s)   & 
     958               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1)      ) & 
    927959               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i & 
    928960               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 

  • NEMO/trunk/src/ICE/iceupdate.F90

    r13643 r14005  
    9494      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
    9595      !!--------------------------------------------------------------------- 
    96       IF( ln_timing )   CALL timing_start('ice_update') 
     96      IF( ln_timing )   CALL timing_start('iceupdate') 
    9797 
    9898      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    154154         ! ice-ocean  mass flux 
    155155         wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    156             &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
     156            &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
    157157          
    158158         ! snw-ocean mass flux 
    … …  
    160160          
    161161         ! total mass flux at the ocean/ice interface 
    162          fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model 
    163          emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux 
     162         fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model 
     163         emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux 
    164164 
    165165         ! Salt flux at the ocean surface       
    … …  
    172172         snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    173173         !                                               ! new mass per unit area 
    174          snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  
     174         snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) )  
    175175         !                                               ! time evolution of snow+ice mass 
    176176         snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice 
    … …  
    286286      IF( ln_icectl         )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    287287      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    288       IF( ln_timing         )   CALL timing_stop   ('ice_update')                                      ! timing 
     288      IF( ln_timing         )   CALL timing_stop   ('iceupdate')                                       ! timing 
    289289      ! 
    290290   END SUBROUTINE ice_update_flx 
    … …  
    324324      REAL(wp) ::   zflagi                          !   -      - 
    325325      !!--------------------------------------------------------------------- 
    326       IF( ln_timing )   CALL timing_start('ice_update_tau') 
     326      IF( ln_timing )   CALL timing_start('ice_update') 
    327327 
    328328      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    376376      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition 
    377377      ! 
    378       IF( ln_timing )   CALL timing_stop('ice_update_tau') 
     378      IF( ln_timing )   CALL timing_stop('ice_update') 
    379379      !   
    380380   END SUBROUTINE ice_update_tau 

  • NEMO/trunk/src/ICE/icevar.F90

    r13472 r14005  
    236236      z1_zhmax =  1._wp / hi_max(jpl)                
    237237      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
    238          h_i  (:,:,jpl) = zhmax 
     238         h_i   (:,:,jpl) = zhmax 
    239239         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    240240         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 
    … …  
    252252      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
    253253      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 ) 
     254         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
    255255      END WHERE 
    256256      ! 
    … …  
    534534         DO_2D( 1, 1, 1, 1 ) 
    535535            ! update exchanges with ocean 
    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 
     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 
    539540            ! 
    540541            a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 
    … …  
    551552            v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    552553            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) 
    553556            ! 
    554557         END_2D 
    … …  
    635638               psv_i  (ji,jj,jl) = 0._wp 
    636639            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 
    637648         END_2D 
    638649         ! 
    … …  
    643654      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp 
    644655      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp 
    645       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    646       WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok 
    647656      ! 
    648657   END SUBROUTINE ice_var_zapneg 
    … …  
    675684      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    676685      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    677       IF( ln_pnd_LEV ) THEN 
     686      IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    678687         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    679688         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 

  • NEMO/trunk/src/ICE/icewri.F90

    r13472 r14005  
    160160      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 
    161161      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories 
     162      IF( iom_use('icevpnd_cat' ) )   CALL iom_put( 'icevpnd_cat' ,   v_ip         * zmsk00l                                   ) ! melt pond volume for categories 
    162163      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 
    163164      IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 
    164       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
     165      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac per ice area for categories 
    165166      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    166167      IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories 
Note: See TracChangeset for help on using the changeset viewer.