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

Changeset 14302


Ignore:
Timestamp:
2021-01-14T17:21:39+01:00 (3 years ago)
Author:
edblockley
Message:

Porting topographic melt-ponds from Martin's branch branches/2020/SI3-05_MeltPonds_topo@14301

Location:
NEMO/branches/UKMO/NEMO_4.0.4_GO8_package
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-ice-BASIC-10D.xml

    r12848 r14302  
    55============================================================================================================ 
    66=                                           output files definition                                        = 
    7 =                                            Define your own files for lim3                                = 
     7=                                      Define your own files for sea ice                                   = 
    88=                                         put the variables you want...                                    = 
    99============================================================================================================ 
    … …  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <!-- sst_m is always the potential temperature even when using teos10 --> 
    3132       <field field_ref="sst_m_pot"        name="sst_m_pot"  /> 
    … …  
    3637       <field field_ref="icettop"          name="sittop" /> 
    3738       <field field_ref="icetbot"          name="sitbot" /> 
     39       <field field_ref="icetsni"          name="sitsni" /> 
     40 
     41       <!-- ponds --> 
     42       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     43       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     44       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     45       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3846        
    3947       <!-- momentum --> 
    … …  
    8593       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8694       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     95       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     96       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8797       <field field_ref="snwtemp_cat"      name="sntemcat"/> 
    8898 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-ice-BASIC-1M.xml

    r12848 r14302  
    55============================================================================================================ 
    66=                                           output files definition                                        = 
    7 =                                            Define your own files for lim3                                = 
     7=                                      Define your own files for sea ice                                   = 
    88=                                         put the variables you want...                                    = 
    99============================================================================================================ 
    … …  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <!-- sst_m is always the potential temperature even when using teos10 --> 
    3132       <field field_ref="sst_m_pot"        name="sst_m_pot"  /> 
    … …  
    3637       <field field_ref="icettop"          name="sittop" /> 
    3738       <field field_ref="icetbot"          name="sitbot" /> 
     39       <field field_ref="icetsni"          name="sitsni" /> 
     40 
     41       <!-- ponds --> 
     42       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     43       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     44       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     45       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3846        
    3947       <!-- momentum --> 
    … …  
    8593       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8694       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     95       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     96       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8797       <field field_ref="snwtemp_cat"      name="sntemcat"/> 
    8898 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-ice-CMIP6-10D.xml

    r12848 r14302  
    55============================================================================================================ 
    66=                                           output files definition                                        = 
    7 =                                            Define your own files for lim3                                = 
     7=                                      Define your own files for sea ice                                   = 
    88=                                         put the variables you want...                                    = 
    99============================================================================================================ 
    … …  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <!-- sst_m is always the potential temperature even when using teos10 --> 
    3132       <field field_ref="sst_m_pot"        name="sst_m_pot"  /> 
    … …  
    3637       <field field_ref="icettop"          name="sittop" /> 
    3738       <field field_ref="icetbot"          name="sitbot" /> 
     39       <field field_ref="icetsni"          name="sitsni" /> 
     40 
     41       <!-- ponds --> 
     42       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     43       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     44       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     45       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3846        
    3947       <!-- momentum --> 
    … …  
    8593       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8694       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     95       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     96       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8797       <field field_ref="snwtemp_cat"      name="sntemcat"/> 
    88         
     98 
    8999       <!-- mass balance --> 
    90100       <field field_ref="dmithd"           name="sidmassth"        /> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/METO_GO/EXPREF/file_def_nemo-ice-CMIP6-1M.xml

    r12848 r14302  
    55============================================================================================================ 
    66=                                           output files definition                                        = 
    7 =                                            Define your own files for lim3                                = 
     7=                                      Define your own files for sea ice                                   = 
    88=                                         put the variables you want...                                    = 
    99============================================================================================================ 
    … …  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <!-- sst_m is always the potential temperature even when using teos10 --> 
    3132       <field field_ref="sst_m_pot"        name="sst_m_pot"  /> 
    … …  
    3637       <field field_ref="icettop"          name="sittop" /> 
    3738       <field field_ref="icetbot"          name="sitbot" /> 
     39       <field field_ref="icetsni"          name="sitsni" /> 
     40 
     41       <!-- ponds --> 
     42       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     43       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     44       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     45       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3846        
    3947       <!-- momentum --> 
    … …  
    8593       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8694       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     95       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     96       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8797       <field field_ref="snwtemp_cat"      name="sntemcat"/> 
    88         
     98 
    8999       <!-- mass balance --> 
    90100       <field field_ref="dmithd"           name="sidmassth"        /> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml

    r9896 r14302  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <field field_ref="sst_m"            name="sst_m"  /> 
    3132       <field field_ref="sss_m"            name="sss_m"  /> 
    … …  
    3637       <field field_ref="icetbot"          name="sitbot" /> 
    3738       <field field_ref="icetsni"          name="sitsni" /> 
     39 
     40       <!-- ponds --> 
     41       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     42       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     43       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     44       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3845        
    3946       <!-- momentum --> 
    … …  
    8592       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8693       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     94       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     95       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8796        
    8897     </file> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/field_def_nemo-ice.xml

    r14078 r14302  
    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      
     54     <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="cm/d" />  
     55     <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="cm/d" />  
     56     <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="cm/d" />  
     57     <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="cm/d" />  
    5358      
    5459     <!-- heat --> 
    … …  
    284289          <field id="snwtemp_cat"  long_name="Snow temperature per category"                     unit="degC"    detect_missing_value="true" /> 
    285290          <field id="icettop_cat"  long_name="Ice/snow surface temperature per category"         unit="degC"    detect_missing_value="true" /> 
    286           <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        />  
     291          <field id="icevpnd_cat"  long_name="Ice melt pond volume per grid area per category"   unit="m"       />  
     292          <field id="iceapnd_cat"  long_name="Ice melt pond grid fraction"                       unit=""        />  
    287293          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
    288294          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    289           <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     295          <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category"           unit=""        />  
    290296          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
    291297          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/cfgs/SHARED/namelist_ice_ref

    r14158 r14302  
    195195!------------------------------------------------------------------------------ 
    196196   ln_pnd            = .true.         !  activate melt ponds or not 
    197       ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 
    198          rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 
    199          rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 
     197      ln_pnd_TOPO    = .true.         !  topographic melt ponds  
     198      ln_pnd_LEV     = .false.        !  level ice melt ponds  
     199         rn_apnd_min =   0.15         !     minimum meltwater fraction contributing to pond growth (TOPO and LEV) 
     200         rn_apnd_max =   0.85         !     maximum meltwater fraction contributing to pond growth (TOPO and LEV) 
    200201         rn_pnd_flush=   0.01         !     pond flushing efficiency (tuning parameter) (LEV) 
    201202      ln_pnd_CST     = .false.        !  constant  melt ponds 
    202203         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC 
    203204         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC 
    204       ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV) 
     205      ln_pnd_lids    = .false.        !  frozen lids on top of the ponds (only for ln_pnd_LEV) 
    205206      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo 
    206207/ 
    … …  
    227228   rn_tms_ini_n     = 270.            !  initial snw temperature     (K), North 
    228229   rn_tms_ini_s     = 270.            !        "            "             South 
    229    rn_apd_ini_n     =   0.2           !  initial pond fraction       (-), North 
     230   rn_apd_ini_n     =   0.0           !  initial pond fraction       (-), North 
    230231   rn_apd_ini_s     =   0.2           !        "            "             South 
    231    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North 
     232   rn_hpd_ini_n     =   0.00          !  initial pond depth          (m), North 
    232233   rn_hpd_ini_s     =   0.05          !        "            "             South 
    233234   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/ice.F90

    r14158 r14302  
    208208   !                                     !!** ice-ponds namelist (namthd_pnd) 
    209209   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    210    LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 
    211    REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds 
    212    REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds 
     210   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 
     211   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Simple melt pond scheme 
     212   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum fraction of melt water contributing to ponds 
     213   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum fraction of melt water contributing to ponds 
     214   REAL(wp), PUBLIC ::   rn_pnd_flush     !: Pond flushing efficiency (tuning parameter) 
    213215   REAL(wp), PUBLIC ::   rn_pnd_flush     !: Pond flushing efficiency (tuning parameter) 
    214216   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    … …  
    309311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K] 
    310312   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 
     313 
     314   ! meltwater arrays to save for melt ponds 
     315   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
     316   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
    311317 
    312318   !!---------------------------------------------------------------------- 
    … …  
    459465      ii = ii + 1 
    460466      ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  & 
     467         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    & 
    461468         &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  & 
    462469         &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  & 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icedyn_adv_pra.F90

    r14075 r14302  
    178178               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    179179            END DO 
    180             IF ( ln_pnd_LEV ) THEN 
     180            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    181181               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
    182182               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
    … …  
    214214            END DO 
    215215            ! 
    216             IF ( ln_pnd_LEV ) THEN 
     216            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    217217               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    218218               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    … …  
    249249                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    250250            END DO 
    251             IF ( ln_pnd_LEV ) THEN 
     251            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    252252               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    253253               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    … …  
    278278         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
    279279            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
    280          IF ( ln_pnd_LEV ) THEN 
     280         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    281281            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
    282282               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
    … …  
    302302               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    303303            END DO 
    304             IF ( ln_pnd_LEV ) THEN 
     304            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    305305               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    306306               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    … …  
    780780                  !                               ! -- check h_ip -- ! 
    781781                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    782                   IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     782                  IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    783783                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    784784                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    … …  
    10271027            END DO 
    10281028            ! 
    1029             IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     1029            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction 
    10301030               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
    10311031                  CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
    … …  
    10691069            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    10701070            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    1071             IF( ln_pnd_LEV ) THEN 
     1071            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    10721072               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
    10731073               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
    … …  
    11371137         END DO 
    11381138         ! 
    1139          IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
     1139         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction 
    11401140            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    11411141            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icedyn_adv_umx.F90

    r14075 r14302  
    341341         ! 
    342342         !== melt ponds ==! 
    343          IF ( ln_pnd_LEV ) THEN 
     343         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    344344            ! concentration 
    345345            zamsk = 1._wp 
    … …  
    361361         ! 
    362362         ! --- Lateral boundary conditions --- ! 
    363          IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     363         IF    ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 
    364364            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 & 
    365365               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
    366          ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     366         ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 
    367367            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 & 
    368368               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
    … …  
    16111611                  !                               ! -- check h_ip -- ! 
    16121612                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1613                   IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1613                  IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    16141614                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    16151615                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icedyn_rdgrft.F90

    r14075 r14302  
    567567               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft  
    568568 
    569                IF ( ln_pnd_LEV ) THEN 
     569               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    570570                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    571571                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    … …  
    604604               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji) 
    605605               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    606                IF ( ln_pnd_LEV ) THEN 
     606               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    607607                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    608608                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    … …  
    701701                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
    702702                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    703                   IF ( ln_pnd_LEV ) THEN 
     703                  IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    704704                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    705705                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/iceitd.F90

    r14075 r14302  
    305305            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    306306               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    307                IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
     307               IF( ln_pnd_LEV .OR. ln_pnd_TOPO )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    308308               h_i_1d(ji) = rn_himin 
    309309            ENDIF 
    … …  
    476476               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    477477               !   
    478                IF ( ln_pnd_LEV ) THEN 
     478               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    479479                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction 
    480480                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icestp.F90

    r14075 r14302  
    414414            wfx_pnd(ji,jj) = 0._wp 
    415415 
     416            dh_i_sum_2d(:,:,:) = 0._wp 
     417            dh_s_mlt_2d(:,:,:) = 0._wp 
     418 
    416419            hfx_thd(ji,jj) = 0._wp   ; 
    417420            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icethd.F90

    r14075 r14302  
    247247            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
    248248                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    249                               CALL ice_thd_pnd                          ! Melt ponds formation 
    250249                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    251250            ENDIF 
    … …  
    268267      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    269268      !                    
     269      IF ( ln_pnd  )          CALL ice_thd_pnd                      ! --- Melt ponds 
     270 
    270271      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
    271272      ! 
    … …  
    536537         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) ) 
    537538         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) ) 
     539         ! ponds 
     540         CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum  (1:npti) , dh_i_sum_2d(:,:,kl) ) 
     541         CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt  (1:npti) , dh_s_mlt_2d(:,:,kl) ) 
    538542         ! SIMIP diagnostics          
    539543         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icethd_pnd.F90

    r14158 r14302  
    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 
     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   REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     44      zhi_min = 0.1_wp        , & ! minimum ice thickness with ponds (m)  
     45      zTd     = 273.0_wp      , & ! temperature difference for freeze-up (C) 
     46      zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     47 
     48   !-------------------------------------------------------------------------- 
     49   !  
     50   ! Pond volume per area budget diags 
     51   !   
     52   ! The idea of diags is the volume of ponds per grid cell area is 
     53   ! 
     54   ! dV/dt = mlt + drn + lid + rnf 
     55   ! mlt   = input from surface melting 
     56   ! drn   = drainage through brine network 
     57   ! lid   = lid growth & melt 
     58   ! rnf   = runoff (water directly removed out of surface melting + overflow) 
     59   ! 
     60   ! In topo mode, the pond water lost because it is in the snow is not included in the budget 
     61   ! 
     62   ! In level mode, all terms are incorporated 
     63 
     64   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 
     65      diag_dvpn_mlt,  &     !! meltwater pond volume input      [m/day] 
     66      diag_dvpn_drn,  &     !! pond volume lost by drainage     [m/day] 
     67      diag_dvpn_lid,  &     !! exchange with lid / refreezing   [m/day] 
     68      diag_dvpn_rnf         !! meltwater pond lost to runoff    [m/day] 
     69       
     70   REAL(wp), ALLOCATABLE, DIMENSION(:) ::   & ! pond volume budget diagnostics (1d) 
     71      diag_dvpn_mlt_1d,  &  !! meltwater pond volume input      [m/day] 
     72      diag_dvpn_drn_1d,  &  !! pond volume lost by drainage     [m/day] 
     73      diag_dvpn_lid_1d,  &  !! exchange with lid / refreezing   [m/day] 
     74      diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day] 
     75   ! 
     76   !-------------------------------------------------------------------------- 
    3977 
    4078   !! * Substitutions 
    … …  
    5290      !!                
    5391      !! ** Purpose :   change melt pond fraction and thickness 
    54       !!                 
     92      !! 
     93      !! Note: Melt ponds affect only radiative transfer for now 
     94      !! 
     95      !!       No heat, no salt. 
     96      !!       The melt water they carry is collected but  
     97      !!       not removed from fw budget or released to the ocean 
     98      !! 
     99      !!       A wfx_pnd has been coded for diagnostic purposes 
     100      !!       It is not fully consistent yet. 
     101      !! 
     102      !!       The current diagnostic lacks a contribution from drainage 
     103      !! 
    55104      !!------------------------------------------------------------------- 
    56       ! 
     105      !! 
     106       
     107      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 
     108      ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 
     109 
     110      diag_dvpn_mlt(:,:) = 0._wp    ;   diag_dvpn_drn(:,:)  = 0._wp 
     111      diag_dvpn_lid(:,:) = 0._wp    ;   diag_dvpn_rnf(:,:)  = 0._wp 
     112      diag_dvpn_mlt_1d(:) = 0._wp   ;   diag_dvpn_drn_1d(:) = 0._wp 
     113      diag_dvpn_lid_1d(:) = 0._wp   ;   diag_dvpn_rnf_1d(:) = 0._wp 
     114 
    57115      SELECT CASE ( nice_pnd ) 
    58116      ! 
    59       CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==! 
     117      CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==! 
    60118         ! 
    61       CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==! 
     119      CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
     120         ! 
     121      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
    62122         ! 
    63123      END SELECT 
    64124      ! 
     125 
     126      IF( iom_use('dvpn_mlt'  ) )   CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 
     127      IF( iom_use('dvpn_lid'  ) )   CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 
     128      IF( iom_use('dvpn_drn'  ) )   CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 
     129      IF( iom_use('dvpn_rnf'  ) )   CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 
     130 
     131      DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 
     132      DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 
     133       
    65134   END SUBROUTINE ice_thd_pnd  
    66135 
    … …  
    75144      !!                to non-zero values when t_su = 0C 
    76145      !! 
    77       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) 
     146      !! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd) 
    78147      !!                 
    79148      !! ** Note   : Coupling with such melt ponds is only radiative 
    … …  
    147216      !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 
    148217      !!  
    149       !! ** Note       :   mostly stolen from CICE 
    150       !! 
    151       !! ** References :   Flocco and Feltham (JGR, 2007) 
    152       !!                   Flocco et al       (JGR, 2010) 
    153       !!                   Holland et al      (J. Clim, 2012) 
     218      !! ** Note       :   Mostly stolen from CICE but not only 
     219      !! 
     220      !! ** References :   Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 
     221      !!                    
    154222      !!------------------------------------------------------------------- 
    155223      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    … …  
    168236      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    169237      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
    170       !! 
    171       INTEGER  ::   ji, jk                            ! loop indices 
     238      REAL(wp) ::   zvold                             ! 
     239      !! 
     240      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     241       
    172242      !!------------------------------------------------------------------- 
    173243      z1_rhow   = 1._wp / rhow  
    174244      z1_aspect = 1._wp / zaspect 
    175245      z1_Tp     = 1._wp / zTp  
    176  
    177       DO ji = 1, npti 
    178          !                                                            !----------------------------------------------------! 
    179          IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
    180             !                                                         !----------------------------------------------------! 
    181             !--- Remove ponds on thin ice or tiny ice fractions 
    182             a_ip_1d(ji)      = 0._wp 
    183             h_ip_1d(ji)      = 0._wp 
    184             h_il_1d(ji)      = 0._wp 
    185             !                                                         !--------------------------------! 
    186          ELSE                                                         ! Case ice thickness >= rn_himin ! 
    187             !                                                         !--------------------------------! 
    188             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
    189             v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
    190             ! 
    191             !------------------! 
    192             ! case ice melting ! 
    193             !------------------! 
    194             ! 
    195             !--- available meltwater for melt ponding ---! 
    196             zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    197             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 
    198             zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
    199             ! 
    200             !--- overflow ---! 
    201             ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    202             !    a_ip_max = zfr_mlt * a_i 
    203             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    204             zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    205             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    206  
    207             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    208             !    h_ip_max = 0.5 * h_i 
    209             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    210             zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    211             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     246       
     247      !----------------------------------------------------------------------------------------------- 
     248      !  Identify grid cells with ice 
     249      !----------------------------------------------------------------------------------------------- 
     250      at_i(:,:) = SUM( a_i, dim=3 ) 
     251      ! 
     252      npti = 0   ;   nptidx(:) = 0 
     253      DO jj = 1, jpj 
     254         DO ji = 1, jpi 
     255            IF ( at_i(ji,jj) > epsi10 ) THEN 
     256               npti = npti + 1 
     257               nptidx( npti ) = (jj - 1) * jpi + ji 
     258            ENDIF 
     259         END DO 
     260      END DO 
     261       
     262      !----------------------------------------------------------------------------------------------- 
     263      ! Prepare 1D arrays 
     264      !----------------------------------------------------------------------------------------------- 
     265 
     266      IF( npti > 0 ) THEN 
     267       
     268         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      ) 
     269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          ) 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          ) 
     271          
     272         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     273         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     274         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     275         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
     276       
     277         DO jl = 1, jpl 
     278          
     279            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d      (1:npti), a_i      (:,:,jl) ) 
     280            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d      (1:npti), h_i      (:,:,jl) ) 
     281            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d     (1:npti), t_su     (:,:,jl) ) 
     282            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     283            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     284            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     285 
     286            CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum    (1:npti), dh_i_sum_2d     (:,:,jl) ) 
     287            CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt    (1:npti), dh_s_mlt_2d     (:,:,jl) ) 
    212288             
    213             !--- Pond growing ---! 
    214             v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    215             ! 
    216             !--- Lid melting ---! 
    217             IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    218             ! 
    219             !--- mass flux ---! 
    220             IF( zdv_mlt > 0._wp ) THEN 
    221                zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    222                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    223                ! 
    224                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    225                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    226                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    227             ENDIF 
    228  
    229             !-------------------! 
    230             ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
    231             !-------------------! 
    232             ! 
    233             zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
    234             ! 
    235             !--- Pond contraction (due to refreezing) ---! 
    236             IF( ln_pnd_lids ) THEN 
    237                ! 
    238                !--- Lid growing and subsequent pond shrinking ---!  
    239                zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    240                   &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    241                 
    242                ! Lid growing 
    243                v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
    244                 
    245                ! Pond shrinking 
    246                v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    247  
    248             ELSE 
    249                ! Pond shrinking 
    250                v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
    251             ENDIF 
    252             ! 
    253             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    254             ! v_ip     = h_ip * a_ip 
    255             ! 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) 
    256             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 
    257             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    258  
    259             !---------------!             
    260             ! Pond flushing ! 
    261             !---------------! 
    262             ! height of top of the pond above sea-level 
    263             zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     289            DO jk = 1, nlay_i 
     290               CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
     291            END DO 
    264292             
    265             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
    266             DO jk = 1, nlay_i 
    267                zsbr = - 1.2_wp                                  & 
    268                   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
    269                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    270                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    271                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
    272             END DO 
    273             zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    274              
    275             ! Do the drainage using Darcy's law 
    276             zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! tunable rn_pnd_flush from hunke et al. (2013) 
    277             zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    278             v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    279              
    280             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    281             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 
    282             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    283  
    284             !--- Corrections and lid thickness ---! 
    285             IF( ln_pnd_lids ) THEN 
    286                !--- retrieve lid thickness from volume ---! 
    287                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    288                ELSE                              ;   h_il_1d(ji) = 0._wp 
    289                ENDIF 
    290                !--- remove ponds if lids are much larger than ponds ---! 
    291                IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     293            !----------------------------------------------------------------------------------------------- 
     294            ! Go for ponds 
     295            !----------------------------------------------------------------------------------------------- 
     296 
     297            DO ji = 1, npti 
     298               !                                                            !----------------------------------------------------! 
     299               IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     300                  !                                                         !----------------------------------------------------! 
     301                  !--- Remove ponds on thin ice or tiny ice fractions 
    292302                  a_ip_1d(ji)      = 0._wp 
    293303                  h_ip_1d(ji)      = 0._wp 
    294304                  h_il_1d(ji)      = 0._wp 
     305                  !                                                         !--------------------------------! 
     306               ELSE                                                         ! Case ice thickness >= rn_himin ! 
     307                  !                                                         !--------------------------------! 
     308                  v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     309                  v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     310                  ! 
     311                  !------------------! 
     312                  ! case ice melting ! 
     313                  !------------------! 
     314                  ! 
     315                  !--- available meltwater for melt ponding ---! 
     316                  zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
     317                  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 
     318                  zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
     319                   
     320                  diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_rdtice                      ! surface melt input diag 
     321                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_rdtice   ! runoff diag 
     322                  ! 
     323                  !--- overflow ---! 
     324                  ! 
     325                  ! 1) area driven overflow 
     326                  ! 
     327                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 
     328                  !    a_ip_max = zfr_mlt * a_i 
     329                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     330                  zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     331                  zvold     = zdv_mlt 
     332                  zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )                      
     333                   
     334                  !  
     335                  ! 2) depth driven overflow 
     336                  ! 
     337                  ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     338                  !    h_ip_max = 0.5 * h_i 
     339                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     340                  zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 
     341                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     342                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_rdtice       ! runoff diag - overflow contribution 
     343             
     344                  !--- Pond growing ---! 
     345                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     346                  ! 
     347                  !--- Lid melting ---! 
     348                  IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     349                  ! 
     350                  !--- mass flux ---! 
     351                  ! MV I would recommend to remove that 
     352                  ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 
     353                   
     354                  IF( zdv_mlt > 0._wp ) THEN 
     355                     zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
     356                     wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     357                     ! 
     358                     zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
     359                     wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     360                     wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     361                  ENDIF 
     362 
     363                  !-------------------! 
     364                  ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     365                  !-------------------! 
     366                  ! 
     367                  zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
     368                  ! 
     369                  !--- Pond contraction (due to refreezing) ---! 
     370                  zvold       = v_ip_1d(ji) ! for diag 
     371 
     372                  IF( ln_pnd_lids ) THEN 
     373                     ! 
     374                     !--- Lid growing and subsequent pond shrinking ---!  
     375                     zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     376                        &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
     377                
     378                     ! Lid growing 
     379                     v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
     380                
     381                     ! Pond shrinking 
     382                     v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
     383                      
     384                  ELSE 
     385                   
     386                     ! Pond shrinking 
     387                     v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
     388                      
     389                  ENDIF 
     390 
     391                  diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_rdtice   ! shrinking counted as lid diagnostic 
     392 
     393                  ! 
     394                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     395                  ! v_ip     = h_ip * a_ip 
     396                  ! 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) 
     397                  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 
     398                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     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) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     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                     ! MV linear expression more consistent & simpler              zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
     414                     ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     415                  END DO 
     416                  zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     417             
     418                  ! Do the drainage using Darcy's law 
     419                  zdv_flush   = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 
     420                  zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
     421                  ! zdv_flush   = 0._wp ! MV remove pond drainage for now 
     422                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     423                   
     424                  diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_rdtice   ! shrinking counted as lid diagnostic 
     425             
     426                  ! MV --- why pond drainage does not give back water into freshwater flux ?             
     427                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     428                  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 
     429                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     430 
     431                  !--- Corrections and lid thickness ---! 
     432                  IF( ln_pnd_lids ) THEN 
     433                     !--- retrieve lid thickness from volume ---! 
     434                     IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     435                     ELSE                              ;   h_il_1d(ji) = 0._wp 
     436                     ENDIF 
     437                     !--- remove ponds if lids are much larger than ponds ---! 
     438                     IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     439                        a_ip_1d(ji)      = 0._wp 
     440                        h_ip_1d(ji)      = 0._wp 
     441                        h_il_1d(ji)      = 0._wp 
     442                     ENDIF 
     443                  ENDIF 
     444                  ! 
    295445               ENDIF 
     446                
     447            END DO ! ji 
     448 
     449            !----------------------------------------------------------------------------------------------- 
     450            ! Retrieve 2D arrays 
     451            !----------------------------------------------------------------------------------------------- 
     452             
     453            v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
     454            v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
     455             
     456            CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     457            CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     458            CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     459            CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d     (1:npti), v_ip     (:,:,jl) ) 
     460            CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d     (1:npti), v_il     (:,:,jl) ) 
     461            DO jk = 1, nlay_i 
     462               CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
     463            END DO 
     464    
     465         END DO ! ji 
     466                   
     467         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 
     468         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        ) 
     469         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
     470          
     471         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     472         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     473         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     474         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                   
     475                            
     476      ! 
     477      ENDIF 
     478       
     479   END SUBROUTINE pnd_LEV 
     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 
     513      ! local variables 
     514      REAL(wp) :: & 
     515         zdHui,   &      ! change in thickness of ice lid (m) 
     516         zomega,  &      ! conduction 
     517         zdTice,  &      ! temperature difference across ice lid (C) 
     518         zdvice,  &      ! change in ice volume (m) 
     519         zTavg,   &      ! mean surface temperature across categories (C) 
     520         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2) 
     521         zTp,     &      ! pond freezing temperature (C) 
     522         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3) 
     523         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding 
     524         z1_rhow, &      ! inverse water density 
     525         zv_pnd  , &     ! volume of meltwater contributing to ponds 
     526         zv_mlt          ! total amount of meltwater produced 
     527 
     528      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp_ini, & !! total melt pond water available before redistribution and drainage 
     529                                        zvolp    , & !! total melt pond water volume 
     530                                        zvolp_res    !! remaining melt pond water available after drainage 
     531                                         
     532      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i 
     533 
     534      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     535 
     536      INTEGER  ::   i_test 
     537 
     538      ! Note 
     539      ! equivalent for CICE translation 
     540      ! a_ip      -> apond 
     541      ! a_ip_frac -> apnd 
     542       
     543      !--------------------------------------------------------------- 
     544      ! Initialise 
     545      !--------------------------------------------------------------- 
     546 
     547      ! Parameters & constants (move to parameters) 
     548      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3) 
     549      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline) 
     550      z1_rhow   = 1._wp / rhow  
     551 
     552      ! Set required ice variables (hard-coded here for now) 
     553!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?)  
     554       
     555      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 
     556      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 
     557      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 
     558      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 
     559       
     560      ! thickness 
     561      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 
     562      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp 
     563      END WHERE 
     564      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 
     565       
     566      !--------------------------------------------------------------- 
     567      ! Change 2D to 1D 
     568      !--------------------------------------------------------------- 
     569      ! MV  
     570      ! a less computing-intensive version would have 2D-1D passage here 
     571      ! use what we have in iceitd.F90 (incremental remapping) 
     572 
     573      !-------------------------------------------------------------- 
     574      ! Collect total available pond water volume 
     575      !-------------------------------------------------------------- 
     576      ! Assuming that meltwater (+rain in principle) runsoff the surface 
     577      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 
     578      ! I cite her words, they are very talkative 
     579      ! "grid cells with very little ice cover (and hence more open water area)  
     580      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 
     581      ! "This results in the same runoff fraction r for each ice category within a grid cell" 
     582       
     583      zvolp(:,:) = 0._wp 
     584 
     585      DO jl = 1, jpl 
     586         DO jj = 1, jpj 
     587            DO ji = 1, jpi 
     588                  
     589               IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
     590             
     591                  !--- Available and contributing meltwater for melt ponding ---! 
     592                  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 
     593                     &    * z1_rhow * a_i(ji,jj,jl) 
     594                      ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 
     595                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds 
     596                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl 
     597 
     598                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice                     ! diags 
     599                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice    
     600 
     601                  !--- Create possible new ponds 
     602                  ! if pond does not exist, create new pond over full ice area 
     603                  !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 
     604                  IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 
     605                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl) 
     606                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE) 
     607                  ENDIF 
     608                   
     609                  !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 
     610                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness 
     611                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     612                   
     613                  !--- Total available pond water volume (pre-existing + newly produced)j 
     614                  zvolp(ji,jj)   = zvolp(ji,jj)   + v_ip(ji,jj,jl)  
     615!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 
     616                    
     617               ENDIF ! a_i 
     618 
     619            END DO! jl             
     620         END DO ! jj 
     621      END DO ! ji 
     622 
     623      zvolp_ini(:,:) = zvolp(:,:) 
     624                   
     625      !-------------------------------------------------------------- 
     626      ! Redistribute and drain water from ponds 
     627      !--------------------------------------------------------------    
     628      CALL ice_thd_pnd_area( zvolp, zvolp_res ) 
     629                                    
     630      !-------------------------------------------------------------- 
     631      ! Melt pond lid growth and melt 
     632      !--------------------------------------------------------------    
     633       
     634      IF( ln_pnd_lids ) THEN 
     635 
     636         DO jj = 1, jpj 
     637            DO ji = 1, jpi 
     638 
     639            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp_ini(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  ! line valid for temp in C  
     691                        zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0    
     692                        zomega = rcnd_i * zdTice / zrhoi_L 
     693                        zdHui  = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 
     694                               - v_il(ji,jj,jl) / a_i(ji,jj,jl) 
     695                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     696                   
     697                        IF ( zdvice > epsi10 ) THEN 
     698                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     699                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     700!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice 
     701!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice 
     702                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     703                            
     704                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag 
     705                            
     706                        ENDIF 
     707                   
     708                     ENDIF ! Tsfcn(i,j,n) 
     709 
     710                  !---------------------------------------------------------------- 
     711                  ! Freeze new lids 
     712                  !---------------------------------------------------------------- 
     713                  !  upper ice layer begins to form 
     714                  ! note: albedo does not change 
     715 
     716                  ELSE ! v_il < epsi10 
     717                     
     718                     ! thickness of newly formed ice 
     719                     ! the surface temperature of a meltpond is the same as that 
     720                     ! of the ice underneath (0C), and the thermodynamic surface  
     721                     ! flux is the same 
     722                      
     723                     !!! we need net surface energy flux, excluding conduction 
     724                     !!! fsurf is summed over categories in CICE 
     725                     !!! we have the category-dependent flux, let us use it ? 
     726                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                      
     727                     zdHui  = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp ) 
     728                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     729 
     730                     IF ( zdvice > epsi10 ) THEN 
     731 
     732                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     733                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     734                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag 
     735 
     736!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice 
     737!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice 
     738 
     739                        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 
     740                     ENDIF 
     741                
     742                  ENDIF  ! v_il 
     743             
     744               END DO ! jl 
     745 
     746            ELSE  ! remove ponds on thin ice 
     747 
     748               v_ip(ji,jj,:) = 0._wp 
     749               v_il(ji,jj,:) = 0._wp 
     750!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 
     751!                 zvolp(ji,jj)    = 0._wp          
     752 
    296753            ENDIF 
    297             ! 
    298          ENDIF 
    299           
    300       END DO 
    301       ! 
    302    END SUBROUTINE pnd_LEV 
    303  
     754 
     755      END DO   ! jj 
     756   END DO   ! ji 
     757 
     758      ENDIF ! ln_pnd_lids 
     759 
     760      !--------------------------------------------------------------- 
     761      ! Clean-up variables (probably duplicates what icevar would do) 
     762      !--------------------------------------------------------------- 
     763      ! MV comment 
     764      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 
     765      ! icevar should recompute all other variables (if needed at all) 
     766 
     767      DO jl = 1, jpl 
     768 
     769         DO jj = 1, jpj 
     770            DO ji = 1, jpi 
     771 
     772!              ! zap lids on small ponds 
     773!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 
     774!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN 
     775!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 
     776!              ENDIF 
     777       
     778               ! recalculate equivalent pond variables 
     779               IF ( a_ip(ji,jj,jl) > epsi10) THEN 
     780                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 
     781                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     782                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     783               ENDIF 
     784!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar 
     785!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar 
     786!              ENDIF 
     787                
     788            END DO   ! ji 
     789         END DO   ! jj 
     790 
     791      END DO   ! jl 
     792 
     793   END SUBROUTINE pnd_TOPO 
     794 
     795    
     796   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 
     797 
     798       !!------------------------------------------------------------------- 
     799       !!                ***  ROUTINE ice_thd_pnd_area *** 
     800       !! 
     801       !! ** Purpose : Given the total volume of available pond water,  
     802       !!              redistribute and drain water 
     803       !! 
     804       !! ** Method 
     805       !! 
     806       !-----------| 
     807       !           | 
     808       !           |-----------| 
     809       !___________|___________|______________________________________sea-level 
     810       !           |           | 
     811       !           |           |---^--------| 
     812       !           |           |   |        | 
     813       !           |           |   |        |-----------|              |------- 
     814       !           |           |   | alfan  |           |              | 
     815       !           |           |   |        |           |--------------| 
     816       !           |           |   |        |           |              | 
     817       !---------------------------v------------------------------------------- 
     818       !           |           |   ^        |           |              | 
     819       !           |           |   |        |           |--------------| 
     820       !           |           |   | betan  |           |              | 
     821       !           |           |   |        |-----------|              |------- 
     822       !           |           |   |        | 
     823       !           |           |---v------- | 
     824       !           |           | 
     825       !           |-----------| 
     826       !           | 
     827       !-----------| 
     828       ! 
     829       !! 
     830       !!------------------------------------------------------------------ 
     831        
     832       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 
     833          zvolp,                                       &  ! total available pond water 
     834          zdvolp                                          ! remaining meltwater after redistribution 
     835 
     836       INTEGER ::  & 
     837          ns,      & 
     838          m_index, & 
     839          permflag 
     840 
     841       REAL (wp), DIMENSION(jpl) :: & 
     842          hicen, & 
     843          hsnon, & 
     844          asnon, & 
     845          alfan, & 
     846          betan, & 
     847          cum_max_vol, & 
     848          reduced_aicen 
     849 
     850       REAL (wp), DIMENSION(0:jpl) :: & 
     851          cum_max_vol_tmp 
     852 
     853       REAL (wp) :: & 
     854          hpond, & 
     855          drain, & 
     856          floe_weight, & 
     857          pressure_head, & 
     858          hsl_rel, & 
     859          deltah, & 
     860          perm, & 
     861          msno 
     862 
     863       REAL (wp), parameter :: & 
     864          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s 
     865 
     866       INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     867 
     868       a_ip(:,:,:) = 0._wp 
     869       h_ip(:,:,:) = 0._wp 
     870  
     871       DO jj = 1, jpj 
     872          DO ji = 1, jpi 
     873  
     874             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     875  
     876        !------------------------------------------------------------------- 
     877        ! initialize 
     878        !------------------------------------------------------------------- 
     879  
     880        DO jl = 1, jpl 
     881  
     882           !---------------------------------------- 
     883           ! compute the effective snow fraction 
     884           !---------------------------------------- 
     885  
     886           IF (a_i(ji,jj,jl) < epsi10)  THEN 
     887              hicen(jl) =  0._wp 
     888              hsnon(jl) =  0._wp 
     889              reduced_aicen(jl) = 0._wp 
     890              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 
     891           ELSE 
     892              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     893              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     894              reduced_aicen(jl) = 1._wp ! n=jpl 
     895  
     896              !js: initial code in NEMO_DEV 
     897              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 
     898              !                     * (-0.024_wp*hicen(jl) + 0.832_wp) 
     899  
     900              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 
     901              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &  
     902                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 
     903  
     904              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical) 
     905              ! MV should check whether this makes sense to have the same effective snow fraction in here 
     906              ! OLI: it probably doesn't 
     907           END IF 
     908  
     909  ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     910  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     911  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     912  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     913  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     914  ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
     915  ! alfan = 60% of the ice volume) in each category lies above the reference line, 
     916  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     917  
     918  ! MV: 
     919  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
     920  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
     921  
     922  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 
     923  
     924           alfan(jl) = 0.6 * hicen(jl) 
     925           betan(jl) = 0.4 * hicen(jl) 
     926  
     927           cum_max_vol(jl)     = 0._wp 
     928           cum_max_vol_tmp(jl) = 0._wp 
     929  
     930        END DO ! jpl 
     931  
     932        cum_max_vol_tmp(0) = 0._wp 
     933        drain = 0._wp 
     934        zdvolp(ji,jj) = 0._wp 
     935  
     936        !---------------------------------------------------------- 
     937        ! Drain overflow water, update pond fraction and volume 
     938        !---------------------------------------------------------- 
     939  
     940        !-------------------------------------------------------------------------- 
     941        ! the maximum amount of water that can be contained up to each ice category 
     942        !-------------------------------------------------------------------------- 
     943        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     944        ! Then the excess volume cum_max_vol(jl) drains out of the system 
     945        ! It should be added to wfx_pnd_out 
     946  
     947        DO jl = 1, jpl-1 ! last category can not hold any volume 
     948  
     949           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 
     950  
     951              ! total volume in level including snow 
     952              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 
     953                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 
     954  
     955              ! subtract snow solid volumes from lower categories in current level 
     956              DO ns = 1, jl 
     957                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 
     958                    - rhos/rhow * &     ! free air fraction that can be filled by water 
     959                      asnon(ns)  * &    ! effective areal fraction of snow in that category 
     960                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 
     961              END DO 
     962  
     963           ELSE ! assume higher categories unoccupied 
     964              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 
     965           END IF 
     966           !IF (cum_max_vol_tmp(jl) < z0) THEN 
     967           !   CALL abort_ice('negative melt pond volume') 
     968           !END IF 
     969        END DO 
     970        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     971        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     972  
     973        !---------------------------------------------------------------- 
     974        ! is there more meltwater than can be held in the floe? 
     975        !---------------------------------------------------------------- 
     976        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 
     977           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 
     978           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 
     979  
     980           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice) 
     981            
     982           zdvolp(ji,jj) = drain         ! this is the drained water 
     983           IF (zvolp(ji,jj) < epsi10) THEN 
     984              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     985              zvolp(ji,jj) = 0._wp 
     986           END IF 
     987        END IF 
     988  
     989        ! height and area corresponding to the remaining volume 
     990        ! routine leaves zvolp unchanged 
     991        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     992  
     993        DO jl = 1, m_index 
     994           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 
     995           !                                         !  volume instead, no ? 
     996           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2 
     997           a_ip(ji,jj,jl) = reduced_aicen(jl) 
     998           ! in practise, pond fraction depends on the empirical snow fraction 
     999           ! so in turn on ice thickness 
     1000        END DO 
     1001        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1002  
     1003        !------------------------------------------------------------------------ 
     1004        ! Drainage through brine network (permeability) 
     1005        !------------------------------------------------------------------------ 
     1006        !!! drainage due to ice permeability - Darcy's law 
     1007  
     1008        ! sea water level 
     1009        msno = 0._wp  
     1010        DO jl = 1 , jpl 
     1011          msno = msno + v_s(ji,jj,jl) * rhos 
     1012        END DO 
     1013        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 
     1014        hsl_rel = floe_weight / rau0 & 
     1015                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 
     1016  
     1017        deltah = hpond - hsl_rel 
     1018        pressure_head = grav * rau0 * max(deltah, 0._wp) 
     1019  
     1020        ! drain if ice is permeable 
     1021        permflag = 0 
     1022  
     1023        IF (pressure_head > 0._wp) THEN 
     1024           DO jl = 1, jpl-1 
     1025              IF ( hicen(jl) /= 0._wp ) THEN 
     1026  
     1027              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2 
     1028  
     1029                 perm = 0._wp ! MV ugly dummy patch 
     1030                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
     1031                 IF (perm > 0._wp) permflag = 1 
     1032  
     1033                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 
     1034                                          (viscosity*hicen(jl)) 
     1035                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 
     1036                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 
     1037  
     1038                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 
     1039                  
     1040                 IF (zvolp(ji,jj) < epsi10) THEN 
     1041                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     1042                    zvolp(ji,jj) = 0._wp 
     1043                 END IF 
     1044             END IF 
     1045          END DO 
     1046  
     1047           ! adjust melt pond dimensions 
     1048           IF (permflag > 0) THEN 
     1049              ! recompute pond depth 
     1050             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     1051              DO jl = 1, m_index 
     1052                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 
     1053                 a_ip(ji,jj,jl) = reduced_aicen(jl) 
     1054              END DO 
     1055              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1056           END IF 
     1057        END IF ! pressure_head 
     1058  
     1059        !------------------------------- 
     1060        ! remove water from the snow 
     1061        !------------------------------- 
     1062        !------------------------------------------------------------------------ 
     1063        ! total melt pond volume in category does not include snow volume 
     1064        ! snow in melt ponds is not melted 
     1065        !------------------------------------------------------------------------ 
     1066         
     1067        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 
     1068        ! how much, so I did not diagnose it 
     1069        ! so if there is a problem here, nobody is going to see it... 
     1070         
     1071  
     1072        ! Calculate pond volume for lower categories 
     1073        DO jl = 1,m_index-1 
     1074           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 
     1075                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
     1076        END DO 
     1077  
     1078        ! Calculate pond volume for highest category = remaining pond volume 
     1079  
     1080        ! The following is completely unclear to Martin at least 
     1081        ! Could we redefine properly and recode in a more readable way ? 
     1082  
     1083        ! m_index = last category with melt pond 
     1084  
     1085        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 
     1086  
     1087        IF (m_index > 1) THEN 
     1088          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 
     1089             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 
     1090          ELSE 
     1091             v_ip(ji,jj,m_index) = 0._wp  
     1092             h_ip(ji,jj,m_index) = 0._wp 
     1093             a_ip(ji,jj,m_index) = 0._wp 
     1094             ! If remaining pond volume is negative reduce pond volume of 
     1095             ! lower category 
     1096             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 
     1097              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) 
     1098          END IF 
     1099        END IF 
     1100  
     1101        DO jl = 1,m_index 
     1102           IF (a_ip(ji,jj,jl) > epsi10) THEN 
     1103               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     1104           ELSE 
     1105              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 
     1106              h_ip(ji,jj,jl) = 0._wp  
     1107              v_ip(ji,jj,jl)  = 0._wp 
     1108              a_ip(ji,jj,jl) = 0._wp 
     1109           END IF 
     1110        END DO 
     1111        DO jl = m_index+1, jpl 
     1112           h_ip(ji,jj,jl) = 0._wp  
     1113           a_ip(ji,jj,jl) = 0._wp  
     1114           v_ip(ji,jj,jl) = 0._wp  
     1115        END DO 
     1116         
     1117           ENDIF 
     1118        END DO ! ji 
     1119     END DO ! jj 
     1120 
     1121    END SUBROUTINE ice_thd_pnd_area 
     1122 
     1123 
     1124    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 
     1125       !!------------------------------------------------------------------- 
     1126       !!                ***  ROUTINE ice_thd_pnd_depth  *** 
     1127       !! 
     1128       !! ** Purpose :   Compute melt pond depth 
     1129       !!------------------------------------------------------------------- 
     1130 
     1131       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     1132          aicen, & 
     1133          asnon, & 
     1134          hsnon, & 
     1135          alfan, & 
     1136          cum_max_vol 
     1137 
     1138       REAL (wp), INTENT(IN) :: & 
     1139          zvolp 
     1140 
     1141       REAL (wp), INTENT(OUT) :: & 
     1142          hpond 
     1143 
     1144       INTEGER, INTENT(OUT) :: & 
     1145          m_index 
     1146 
     1147       INTEGER :: n, ns 
     1148 
     1149       REAL (wp), DIMENSION(0:jpl+1) :: & 
     1150          hitl, & 
     1151          aicetl 
     1152 
     1153       REAL (wp) :: & 
     1154          rem_vol, & 
     1155          area, & 
     1156          vol, & 
     1157          tmp, & 
     1158          z0   = 0.0_wp 
     1159 
     1160       !---------------------------------------------------------------- 
     1161       ! hpond is zero if zvolp is zero - have we fully drained? 
     1162       !---------------------------------------------------------------- 
     1163 
     1164       IF (zvolp < epsi10) THEN 
     1165        hpond = z0 
     1166        m_index = 0 
     1167       ELSE 
     1168 
     1169        !---------------------------------------------------------------- 
     1170        ! Calculate the category where water fills up to 
     1171        !---------------------------------------------------------------- 
     1172 
     1173        !----------| 
     1174        !          | 
     1175        !          | 
     1176        !          |----------|                                     -- -- 
     1177        !__________|__________|_________________________________________ ^ 
     1178        !          |          |             rem_vol     ^                | Semi-filled 
     1179        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     1180        !          |          |          |              | 
     1181        !          |          |          |              |hpond 
     1182        !          |          |          |----------|   |     |------- 
     1183        !          |          |          |          |   |     | 
     1184        !          |          |          |          |---v-----| 
     1185        !          |          | m_index  |          |         | 
     1186        !------------------------------------------------------------- 
     1187 
     1188        m_index = 0  ! 1:m_index categories have water in them 
     1189        DO n = 1, jpl 
     1190           IF (zvolp <= cum_max_vol(n)) THEN 
     1191              m_index = n 
     1192              IF (n == 1) THEN 
     1193                 rem_vol = zvolp 
     1194              ELSE 
     1195                 rem_vol = zvolp - cum_max_vol(n-1) 
     1196              END IF 
     1197              exit ! to break out of the loop 
     1198           END IF 
     1199        END DO 
     1200        m_index = min(jpl-1, m_index) 
     1201 
     1202        !---------------------------------------------------------------- 
     1203        ! semi-filled layer may have m_index different snow in it 
     1204        !---------------------------------------------------------------- 
     1205 
     1206        !-----------------------------------------------------------  ^ 
     1207        !                                                             |  alfan(m_index+1) 
     1208        !                                                             | 
     1209        !hitl(3)-->                             |----------|          | 
     1210        !hitl(2)-->                |------------| * * * * *|          | 
     1211        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     1212        !hitl(0)-->-------------------------------------------------  |  ^ 
     1213        !                various snow from lower categories          |  |alfa(m_index) 
     1214 
     1215        ! hitl - heights of the snow layers from thinner and current categories 
     1216        ! aicetl - area of each snow depth in this layer 
     1217 
     1218        hitl(:) = z0 
     1219        aicetl(:) = z0 
     1220        DO n = 1, m_index 
     1221           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     1222                                  alfan(m_index+1) - alfan(m_index)), z0) 
     1223           aicetl(n) = asnon(n) 
     1224 
     1225           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     1226        END DO 
     1227 
     1228        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     1229        aicetl(m_index+1) = z0 
     1230 
     1231        !---------------------------------------------------------------- 
     1232        ! reorder array according to hitl 
     1233        ! snow heights not necessarily in height order 
     1234        !---------------------------------------------------------------- 
     1235 
     1236        DO ns = 1, m_index+1 
     1237           DO n = 0, m_index - ns + 1 
     1238              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     1239                 tmp = hitl(n) 
     1240                 hitl(n) = hitl(n+1) 
     1241                 hitl(n+1) = tmp 
     1242                 tmp = aicetl(n) 
     1243                 aicetl(n) = aicetl(n+1) 
     1244                 aicetl(n+1) = tmp 
     1245              END IF 
     1246           END DO 
     1247        END DO 
     1248 
     1249        !---------------------------------------------------------------- 
     1250        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     1251        !---------------------------------------------------------------- 
     1252 
     1253        !hitl(3)---------------------------------------------------------------- 
     1254        !                                                       | * * * * * * * * 
     1255        !                                                       |* * * * * * * * * 
     1256        !hitl(2)---------------------------------------------------------------- 
     1257        !                                    | * * * * * * * *  | * * * * * * * * 
     1258        !                                    |* * * * * * * * * |* * * * * * * * * 
     1259        !hitl(1)---------------------------------------------------------------- 
     1260        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
     1261        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
     1262        !hitl(0)---------------------------------------------------------------- 
     1263        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3) 
     1264 
     1265        ! move up over layers incrementing volume 
     1266        DO n = 1, m_index+1 
     1267 
     1268           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     1269                (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     1270 
     1271           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     1272 
     1273           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     1274              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 
     1275 
     1276              exit 
     1277           ELSE  ! still in sub-layer below the sub-layer with the depth 
     1278              rem_vol = rem_vol - vol 
     1279           END IF 
     1280 
     1281        END DO 
     1282 
     1283       END IF 
     1284 
     1285    END SUBROUTINE ice_thd_pnd_depth 
     1286 
     1287 
     1288    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 
     1289       !!------------------------------------------------------------------- 
     1290       !!                ***  ROUTINE ice_thd_pnd_perm *** 
     1291       !! 
     1292       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     1293       !!                and its permeability 
     1294       !!------------------------------------------------------------------- 
     1295 
     1296       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     1297          ticen,  &     ! internal ice temperature (K) 
     1298          salin         ! salinity (ppt)     !js: ppt according to cice 
     1299 
     1300       REAL (wp), INTENT(OUT) :: & 
     1301          perm      ! permeability 
     1302 
     1303       REAL (wp) ::   & 
     1304          Sbr       ! brine salinity 
     1305 
     1306       REAL (wp), DIMENSION(nlay_i) ::   & 
     1307          Tin, &    ! ice temperature 
     1308          phi       ! liquid fraction 
     1309 
     1310       INTEGER :: k 
     1311 
     1312       !----------------------------------------------------------------- 
     1313       ! Compute ice temperatures from enthalpies using quadratic formula 
     1314       !----------------------------------------------------------------- 
     1315 
     1316       DO k = 1,nlay_i 
     1317          Tin(k) = ticen(k) - rt0   !js: from K to degC 
     1318       END DO 
     1319 
     1320       !----------------------------------------------------------------- 
     1321       ! brine salinity and liquid fraction 
     1322       !----------------------------------------------------------------- 
     1323 
     1324       DO k = 1, nlay_i 
     1325        
     1326          ! Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 
     1327          ! Best expression to date is that one 
     1328          Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 
     1329          phi(k) = salin(k) / Sbr 
     1330           
     1331       END DO 
     1332 
     1333       !----------------------------------------------------------------- 
     1334       ! permeability 
     1335       !----------------------------------------------------------------- 
     1336 
     1337       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 
     1338 
     1339   END SUBROUTINE ice_thd_pnd_perm 
     1340    
     1341    
     1342!---------------------------------------------------------------------------------------------------------------------- 
    3041343 
    3051344   SUBROUTINE ice_thd_pnd_init  
    … …  
    3181357      !! 
    3191358      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 
    320          &                          ln_pnd_CST , rn_apnd, rn_hpnd,         & 
     1359         &                          ln_pnd_CST , rn_apnd, rn_hpnd,                       & 
     1360         &                          ln_pnd_TOPO ,                                        & 
    3211361         &                          ln_pnd_lids, ln_pnd_alb 
    3221362      !!------------------------------------------------------------------- 
    … …  
    3361376         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    3371377         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd 
     1378         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO 
    3381379         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV 
    3391380         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min 
    … …  
    3521393      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    3531394      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF 
     1395      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF 
    3541396      IF( ioptio /= 1 )   & 
    355          & 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)' ) 
     1397         & 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)' ) 
    3561398      ! 
    3571399      SELECT CASE( nice_pnd ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icevar.F90

    r14075 r14302  
    565565            END DO 
    566566         END DO 
     567 
     568         !----------------------------------------------------------------- 
     569         ! zap small ponds 
     570         !----------------------------------------------------------------- 
     571         DO jj = 1, jpj 
     572            DO ji = 1, jpi 
     573               IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN 
     574                  a_ip(ji,jj,jl)      = 0._wp 
     575                  a_ip_frac(ji,jj,jl) = 0._wp 
     576                  v_ip(ji,jj,jl)      = 0._wp 
     577                  h_ip(ji,jj,jl)      = 0._wp 
     578                  v_il(ji,jj,jl)      = 0._wp 
     579                  h_il(ji,jj,jl)      = 0._wp 
     580                ENDIF 
     581            END DO 
     582         END DO 
    567583         ! 
    568584      END DO  
    … …  
    696712      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    697713      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    698       IF( ln_pnd_LEV ) THEN 
     714      IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    699715         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    700716         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/ICE/icewri.F90

    r14078 r14302  
    164164      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l                                   ) ! brine volume 
    165165      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories 
     166      IF( iom_use('icevpnd_cat' ) )   CALL iom_put( 'icevpnd_cat' ,   v_ip         * zmsk00l                                   ) ! melt pond volume for categories 
    166167      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 
    167168      IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 
    168       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
     169      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac per ice area for categories 
    169170      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    170171      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.