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 14671 for NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_minor/src/ICE – NEMO

Ignore:
Timestamp:
2021-04-01T13:34:55+02:00 (3 years ago)
Author:
dancopsey
Message:

Merged in up to revision 14474 of the GO8_package branch

Location:
NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_minor/src/ICE
Files:
13 edited

Legend:

Unmodified
Added
Removed

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

    r14075 r14671  
    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   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    214216   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
    … …  
    308310   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K] 
    309311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 
     312 
     313   ! meltwater arrays to save for melt ponds 
     314   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
     315   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
    310316 
    311317   !!---------------------------------------------------------------------- 
    … …  
    458464      ii = ii + 1 
    459465      ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  & 
     466         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    & 
    460467         &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  & 
    461468         &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  & 

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

    r14075 r14671  
    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_fix_topo_minor/src/ICE/icedyn_adv_umx.F90

    r14075 r14671  
    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_fix_topo_minor/src/ICE/icedyn_rdgrft.F90

    r14075 r14671  
    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_fix_topo_minor/src/ICE/iceitd.F90

    r14075 r14671  
    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_fix_topo_minor/src/ICE/icerst.F90

    r14075 r14671  
    2727   USE in_out_manager ! I/O manager 
    2828   USE iom            ! I/O manager library 
     29   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2930   USE lib_mpp        ! MPP library 
    3031   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    … …  
    5253      INTEGER, INTENT(in) ::   kt       ! number of iteration 
    5354      ! 
     55      INTEGER             ::   iyear, imonth, iday 
     56      REAL (wp)           ::   zsec 
     57      REAL (wp)           ::   zfjulday 
    5458      CHARACTER(len=20)   ::   clkt     ! ocean time-step define as a character 
    5559      CHARACTER(len=50)   ::   clname   ! ice output restart file name 
    … …  
    6771         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
    6872            ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    69             IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    70             ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     73            IF ( ln_rstdate ) THEN 
     74               zfjulday = fjulday + (2*nn_fsbc+1)*rdt / rday 
     75               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     76               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     77               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     78            ELSE 
     79               IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     80               ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     81               ENDIF 
    7182            ENDIF 
    7283            ! create the file 

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

    r14075 r14671  
    211211      ! --- Ocean time step --- ! 
    212212      !-------------------------! 
    213       IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses 
     213      CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses 
    214214!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    215215      ! 
    … …  
    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_fix_topo_minor/src/ICE/icetab.F90

    r14075 r14671  
    4040      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index 
    4141      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field 
    42       REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field 
     42      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab1d    ! output 1D field 
    4343      ! 
    4444      INTEGER ::   jl, jn, jid, jjd 
    … …  
    6161      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind  ! input index 
    6262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   tab2d    ! input 2D field 
    63       REAL(wp), DIMENSION(ndim1d) , INTENT(  out) ::   tab1d    ! output 1D field 
     63      REAL(wp), DIMENSION(ndim1d) , INTENT(inout) ::   tab1d    ! output 1D field 
    6464      ! 
    6565      INTEGER ::   jn , jid, jjd 
    … …  
    8080      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index 
    8181      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field 
    82       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field 
     82      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab2d     ! output 2D field 
    8383      ! 
    8484      INTEGER ::   jl, jn, jid, jjd 
    … …  
    101101      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind   ! input index 
    102102      REAL(wp), DIMENSION(ndim1d) , INTENT(in   ) ::   tab1d     ! input 1D field 
    103       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   tab2d     ! output 2D field 
     103      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   tab2d     ! output 2D field 
    104104      ! 
    105105      INTEGER ::   jn , jid, jjd 

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

    r14075 r14671  
    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_fix_topo_minor/src/ICE/icethd_pnd.F90

    r14075 r14671  
    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 
    … …  
    145214      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect 
    146215      !! 
    147       !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 
     216      !! ** 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) 
    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!! EWB: bug here - should be divided by at_i(ji,jj) but zTavg is not used anywhere 
     650               zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 
     651          
     652               DO jl = 1, jpl-1 
     653             
     654                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN 
     655                
     656                     !---------------------------------------------------------------- 
     657                     ! Lid melting: floating upper ice layer melts in whole or part 
     658                     !---------------------------------------------------------------- 
     659                     ! Use Tsfc for each category 
     660                     IF ( t_su(ji,jj,jl) > zTp ) THEN 
     661 
     662                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 
     663                         
     664                        IF ( zdvice > epsi10 ) THEN 
     665                         
     666                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice 
     667                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -  
     668                                                                        ! as it is already counted in surface melt 
     669!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar) 
     670!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw) 
     671                      
     672                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 
     673                           ! ice lid melted and category is pond covered 
     674                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl)  
     675!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl)  
     676                              v_il(ji,jj,jl)   = 0._wp 
     677                           ENDIF 
     678                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 
     679                            
     680                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag 
     681                            
     682                        ENDIF 
     683                         
     684                     !---------------------------------------------------------------- 
     685                     ! Freeze pre-existing lid  
     686                     !---------------------------------------------------------------- 
     687 
     688                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 
     689 
     690                        ! differential growth of base of surface floating ice layer 
     691                        ! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0  ! line valid for temp in C  
     692                        zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0    
     693                        zomega = rcnd_i * zdTice / zrhoi_L 
     694                        zdHui  = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 
     695                               - v_il(ji,jj,jl) / a_i(ji,jj,jl) 
     696                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     697                   
     698                        IF ( zdvice > epsi10 ) THEN 
     699                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     700                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     701!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice 
     702!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice 
     703                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     704                            
     705                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag 
     706                            
     707                        ENDIF 
     708                   
     709                     ENDIF ! Tsfcn(i,j,n) 
     710 
     711                  !---------------------------------------------------------------- 
     712                  ! Freeze new lids 
     713                  !---------------------------------------------------------------- 
     714                  !  upper ice layer begins to form 
     715                  ! note: albedo does not change 
     716 
     717                  ELSE ! v_il < epsi10 
     718                     
     719                     ! thickness of newly formed ice 
     720                     ! the surface temperature of a meltpond is the same as that 
     721                     ! of the ice underneath (0C), and the thermodynamic surface  
     722                     ! flux is the same 
     723                      
     724                     !!! we need net surface energy flux, excluding conduction 
     725                     !!! fsurf is summed over categories in CICE 
     726                     !!! we have the category-dependent flux, let us use it ? 
     727                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                      
     728                     zdHui  = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp ) 
     729                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     730 
     731                     IF ( zdvice > epsi10 ) THEN 
     732 
     733                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     734                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     735                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag 
     736 
     737!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice 
     738!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice 
     739 
     740                        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 
     741                     ENDIF 
     742                
     743                  ENDIF  ! v_il 
     744             
     745               END DO ! jl 
     746 
     747            ELSE  ! remove ponds on thin ice 
     748 
     749               v_ip(ji,jj,:) = 0._wp 
     750               v_il(ji,jj,:) = 0._wp 
     751!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 
     752!                 zvolp(ji,jj)    = 0._wp          
     753 
    296754            ENDIF 
    297             ! 
    298          ENDIF 
    299           
    300       END DO 
    301       ! 
    302    END SUBROUTINE pnd_LEV 
    303  
     755 
     756      END DO   ! jj 
     757   END DO   ! ji 
     758 
     759      ENDIF ! ln_pnd_lids 
     760 
     761      !--------------------------------------------------------------- 
     762      ! Clean-up variables (probably duplicates what icevar would do) 
     763      !--------------------------------------------------------------- 
     764      ! MV comment 
     765      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 
     766      ! icevar should recompute all other variables (if needed at all) 
     767 
     768      DO jl = 1, jpl 
     769 
     770         DO jj = 1, jpj 
     771            DO ji = 1, jpi 
     772 
     773!              ! zap lids on small ponds 
     774!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 
     775!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN 
     776!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 
     777!              ENDIF 
     778       
     779               ! recalculate equivalent pond variables 
     780               IF ( a_ip(ji,jj,jl) > epsi10) THEN 
     781                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     782                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     783                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     784               ENDIF 
     785!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar 
     786!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar 
     787!              ENDIF 
     788                
     789            END DO   ! ji 
     790         END DO   ! jj 
     791 
     792      END DO   ! jl 
     793 
     794   END SUBROUTINE pnd_TOPO 
     795 
     796    
     797   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 
     798 
     799       !!------------------------------------------------------------------- 
     800       !!                ***  ROUTINE ice_thd_pnd_area *** 
     801       !! 
     802       !! ** Purpose : Given the total volume of available pond water,  
     803       !!              redistribute and drain water 
     804       !! 
     805       !! ** Method 
     806       !! 
     807       !-----------| 
     808       !           | 
     809       !           |-----------| 
     810       !___________|___________|______________________________________sea-level 
     811       !           |           | 
     812       !           |           |---^--------| 
     813       !           |           |   |        | 
     814       !           |           |   |        |-----------|              |------- 
     815       !           |           |   | alfan  |           |              | 
     816       !           |           |   |        |           |--------------| 
     817       !           |           |   |        |           |              | 
     818       !---------------------------v------------------------------------------- 
     819       !           |           |   ^        |           |              | 
     820       !           |           |   |        |           |--------------| 
     821       !           |           |   | betan  |           |              | 
     822       !           |           |   |        |-----------|              |------- 
     823       !           |           |   |        | 
     824       !           |           |---v------- | 
     825       !           |           | 
     826       !           |-----------| 
     827       !           | 
     828       !-----------| 
     829       ! 
     830       !! 
     831       !!------------------------------------------------------------------ 
     832        
     833       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 
     834          zvolp,                                       &  ! total available pond water 
     835          zdvolp                                          ! remaining meltwater after redistribution 
     836 
     837       INTEGER ::  & 
     838          ns,      & 
     839          m_index, & 
     840          permflag 
     841 
     842       REAL (wp), DIMENSION(jpl) :: & 
     843          hicen, & 
     844          hsnon, & 
     845          asnon, & 
     846          alfan, & 
     847          betan, & 
     848          cum_max_vol, & 
     849          reduced_aicen 
     850 
     851       REAL (wp), DIMENSION(0:jpl) :: & 
     852          cum_max_vol_tmp 
     853 
     854       REAL (wp) :: & 
     855          hpond, & 
     856          drain, & 
     857          floe_weight, & 
     858          pressure_head, & 
     859          hsl_rel, & 
     860          deltah, & 
     861          perm, & 
     862          msno 
     863 
     864       REAL (wp), parameter :: & 
     865          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s 
     866 
     867       INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     868 
     869       a_ip(:,:,:) = 0._wp 
     870       h_ip(:,:,:) = 0._wp 
     871  
     872       DO jj = 1, jpj 
     873          DO ji = 1, jpi 
     874  
     875             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 
     876  
     877        !------------------------------------------------------------------- 
     878        ! initialize 
     879        !------------------------------------------------------------------- 
     880  
     881        DO jl = 1, jpl 
     882  
     883           !---------------------------------------- 
     884           ! compute the effective snow fraction 
     885           !---------------------------------------- 
     886  
     887           IF (a_i(ji,jj,jl) < epsi10)  THEN 
     888              hicen(jl) =  0._wp 
     889              hsnon(jl) =  0._wp 
     890              reduced_aicen(jl) = 0._wp 
     891              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 
     892           ELSE 
     893              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     894              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     895              reduced_aicen(jl) = 1._wp ! n=jpl 
     896  
     897              !js: initial code in NEMO_DEV 
     898              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 
     899              !                     * (-0.024_wp*hicen(jl) + 0.832_wp) 
     900  
     901              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 
     902              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &  
     903                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 
     904  
     905              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical) 
     906              ! MV should check whether this makes sense to have the same effective snow fraction in here 
     907              ! OLI: it probably doesn't 
     908           END IF 
     909  
     910  ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     911  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     912  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     913  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     914  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     915  ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
     916  ! alfan = 60% of the ice volume) in each category lies above the reference line, 
     917  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     918  
     919  ! MV: 
     920  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
     921  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
     922  
     923  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 
     924  
     925           alfan(jl) = 0.6 * hicen(jl) 
     926           betan(jl) = 0.4 * hicen(jl) 
     927  
     928           cum_max_vol(jl)     = 0._wp 
     929           cum_max_vol_tmp(jl) = 0._wp 
     930  
     931        END DO ! jpl 
     932  
     933        cum_max_vol_tmp(0) = 0._wp 
     934        drain = 0._wp 
     935        zdvolp(ji,jj) = 0._wp 
     936  
     937        !---------------------------------------------------------- 
     938        ! Drain overflow water, update pond fraction and volume 
     939        !---------------------------------------------------------- 
     940  
     941        !-------------------------------------------------------------------------- 
     942        ! the maximum amount of water that can be contained up to each ice category 
     943        !-------------------------------------------------------------------------- 
     944        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     945        ! Then the excess volume cum_max_vol(jl) drains out of the system 
     946        ! It should be added to wfx_pnd_out 
     947  
     948        DO jl = 1, jpl-1 ! last category can not hold any volume 
     949  
     950           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 
     951  
     952              ! total volume in level including snow 
     953              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 
     954                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 
     955  
     956              ! subtract snow solid volumes from lower categories in current level 
     957              DO ns = 1, jl 
     958                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 
     959                    - rhos/rhow * &     ! free air fraction that can be filled by water 
     960                      asnon(ns)  * &    ! effective areal fraction of snow in that category 
     961                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 
     962              END DO 
     963  
     964           ELSE ! assume higher categories unoccupied 
     965              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 
     966           END IF 
     967           !IF (cum_max_vol_tmp(jl) < z0) THEN 
     968           !   CALL abort_ice('negative melt pond volume') 
     969           !END IF 
     970        END DO 
     971        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     972        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     973  
     974        !---------------------------------------------------------------- 
     975        ! is there more meltwater than can be held in the floe? 
     976        !---------------------------------------------------------------- 
     977        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 
     978           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 
     979           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 
     980  
     981           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice) 
     982            
     983           zdvolp(ji,jj) = drain         ! this is the drained water 
     984           IF (zvolp(ji,jj) < epsi10) THEN 
     985              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     986              zvolp(ji,jj) = 0._wp 
     987           END IF 
     988        END IF 
     989  
     990        ! height and area corresponding to the remaining volume 
     991        ! routine leaves zvolp unchanged 
     992        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     993  
     994        DO jl = 1, m_index 
     995           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 
     996           !                                         !  volume instead, no ? 
     997           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2 
     998           a_ip(ji,jj,jl) = reduced_aicen(jl) 
     999           ! in practise, pond fraction depends on the empirical snow fraction 
     1000           ! so in turn on ice thickness 
     1001        END DO 
     1002        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1003  
     1004        !------------------------------------------------------------------------ 
     1005        ! Drainage through brine network (permeability) 
     1006        !------------------------------------------------------------------------ 
     1007        !!! drainage due to ice permeability - Darcy's law 
     1008  
     1009        ! sea water level 
     1010        msno = 0._wp  
     1011        DO jl = 1 , jpl 
     1012          msno = msno + v_s(ji,jj,jl) * rhos 
     1013        END DO 
     1014        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 
     1015        hsl_rel = floe_weight / rau0 & 
     1016                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 
     1017  
     1018        deltah = hpond - hsl_rel 
     1019        pressure_head = grav * rau0 * max(deltah, 0._wp) 
     1020  
     1021        ! drain if ice is permeable 
     1022        permflag = 0 
     1023  
     1024        IF (pressure_head > 0._wp) THEN 
     1025           DO jl = 1, jpl-1 
     1026              IF ( hicen(jl) /= 0._wp ) THEN 
     1027  
     1028              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2 
     1029  
     1030                 perm = 0._wp ! MV ugly dummy patch 
     1031                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
     1032                 IF (perm > 0._wp) permflag = 1 
     1033  
     1034                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 
     1035                                          (viscosity*hicen(jl)) 
     1036                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 
     1037                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 
     1038  
     1039                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 
     1040                  
     1041                 IF (zvolp(ji,jj) < epsi10) THEN 
     1042                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     1043                    zvolp(ji,jj) = 0._wp 
     1044                 END IF 
     1045             END IF 
     1046          END DO 
     1047  
     1048           ! adjust melt pond dimensions 
     1049           IF (permflag > 0) THEN 
     1050              ! recompute pond depth 
     1051             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     1052              DO jl = 1, m_index 
     1053                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 
     1054                 a_ip(ji,jj,jl) = reduced_aicen(jl) 
     1055              END DO 
     1056              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1057           END IF 
     1058        END IF ! pressure_head 
     1059  
     1060        !------------------------------- 
     1061        ! remove water from the snow 
     1062        !------------------------------- 
     1063        !------------------------------------------------------------------------ 
     1064        ! total melt pond volume in category does not include snow volume 
     1065        ! snow in melt ponds is not melted 
     1066        !------------------------------------------------------------------------ 
     1067         
     1068        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 
     1069        ! how much, so I did not diagnose it 
     1070        ! so if there is a problem here, nobody is going to see it... 
     1071         
     1072  
     1073        ! Calculate pond volume for lower categories 
     1074        DO jl = 1,m_index-1 
     1075           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 
     1076                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
     1077        END DO 
     1078  
     1079        ! Calculate pond volume for highest category = remaining pond volume 
     1080  
     1081        ! The following is completely unclear to Martin at least 
     1082        ! Could we redefine properly and recode in a more readable way ? 
     1083  
     1084        ! m_index = last category with melt pond 
     1085  
     1086        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 
     1087  
     1088        IF (m_index > 1) THEN 
     1089          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 
     1090             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 
     1091          ELSE 
     1092             v_ip(ji,jj,m_index) = 0._wp  
     1093             h_ip(ji,jj,m_index) = 0._wp 
     1094             a_ip(ji,jj,m_index) = 0._wp 
     1095             ! If remaining pond volume is negative reduce pond volume of 
     1096             ! lower category 
     1097             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 
     1098              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) 
     1099          END IF 
     1100        END IF 
     1101  
     1102        DO jl = 1,m_index 
     1103           IF (a_ip(ji,jj,jl) > epsi10) THEN 
     1104               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     1105           ELSE 
     1106              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 
     1107              h_ip(ji,jj,jl) = 0._wp  
     1108              v_ip(ji,jj,jl)  = 0._wp 
     1109              a_ip(ji,jj,jl) = 0._wp 
     1110           END IF 
     1111        END DO 
     1112        DO jl = m_index+1, jpl 
     1113           h_ip(ji,jj,jl) = 0._wp  
     1114           a_ip(ji,jj,jl) = 0._wp  
     1115           v_ip(ji,jj,jl) = 0._wp  
     1116        END DO 
     1117         
     1118           ENDIF 
     1119        END DO ! ji 
     1120     END DO ! jj 
     1121 
     1122    END SUBROUTINE ice_thd_pnd_area 
     1123 
     1124 
     1125    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 
     1126       !!------------------------------------------------------------------- 
     1127       !!                ***  ROUTINE ice_thd_pnd_depth  *** 
     1128       !! 
     1129       !! ** Purpose :   Compute melt pond depth 
     1130       !!------------------------------------------------------------------- 
     1131 
     1132       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     1133          aicen, & 
     1134          asnon, & 
     1135          hsnon, & 
     1136          alfan, & 
     1137          cum_max_vol 
     1138 
     1139       REAL (wp), INTENT(IN) :: & 
     1140          zvolp 
     1141 
     1142       REAL (wp), INTENT(OUT) :: & 
     1143          hpond 
     1144 
     1145       INTEGER, INTENT(OUT) :: & 
     1146          m_index 
     1147 
     1148       INTEGER :: n, ns 
     1149 
     1150       REAL (wp), DIMENSION(0:jpl+1) :: & 
     1151          hitl, & 
     1152          aicetl 
     1153 
     1154       REAL (wp) :: & 
     1155          rem_vol, & 
     1156          area, & 
     1157          vol, & 
     1158          tmp, & 
     1159          z0   = 0.0_wp 
     1160 
     1161       !---------------------------------------------------------------- 
     1162       ! hpond is zero if zvolp is zero - have we fully drained? 
     1163       !---------------------------------------------------------------- 
     1164 
     1165       IF (zvolp < epsi10) THEN 
     1166        hpond = z0 
     1167        m_index = 0 
     1168       ELSE 
     1169 
     1170        !---------------------------------------------------------------- 
     1171        ! Calculate the category where water fills up to 
     1172        !---------------------------------------------------------------- 
     1173 
     1174        !----------| 
     1175        !          | 
     1176        !          | 
     1177        !          |----------|                                     -- -- 
     1178        !__________|__________|_________________________________________ ^ 
     1179        !          |          |             rem_vol     ^                | Semi-filled 
     1180        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     1181        !          |          |          |              | 
     1182        !          |          |          |              |hpond 
     1183        !          |          |          |----------|   |     |------- 
     1184        !          |          |          |          |   |     | 
     1185        !          |          |          |          |---v-----| 
     1186        !          |          | m_index  |          |         | 
     1187        !------------------------------------------------------------- 
     1188 
     1189        m_index = 0  ! 1:m_index categories have water in them 
     1190        DO n = 1, jpl 
     1191           IF (zvolp <= cum_max_vol(n)) THEN 
     1192              m_index = n 
     1193              IF (n == 1) THEN 
     1194                 rem_vol = zvolp 
     1195              ELSE 
     1196                 rem_vol = zvolp - cum_max_vol(n-1) 
     1197              END IF 
     1198              exit ! to break out of the loop 
     1199           END IF 
     1200        END DO 
     1201        m_index = min(jpl-1, m_index) 
     1202 
     1203        !---------------------------------------------------------------- 
     1204        ! semi-filled layer may have m_index different snow in it 
     1205        !---------------------------------------------------------------- 
     1206 
     1207        !-----------------------------------------------------------  ^ 
     1208        !                                                             |  alfan(m_index+1) 
     1209        !                                                             | 
     1210        !hitl(3)-->                             |----------|          | 
     1211        !hitl(2)-->                |------------| * * * * *|          | 
     1212        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     1213        !hitl(0)-->-------------------------------------------------  |  ^ 
     1214        !                various snow from lower categories          |  |alfa(m_index) 
     1215 
     1216        ! hitl - heights of the snow layers from thinner and current categories 
     1217        ! aicetl - area of each snow depth in this layer 
     1218 
     1219        hitl(:) = z0 
     1220        aicetl(:) = z0 
     1221        DO n = 1, m_index 
     1222           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     1223                                  alfan(m_index+1) - alfan(m_index)), z0) 
     1224           aicetl(n) = asnon(n) 
     1225 
     1226           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     1227        END DO 
     1228 
     1229        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     1230        aicetl(m_index+1) = z0 
     1231 
     1232        !---------------------------------------------------------------- 
     1233        ! reorder array according to hitl 
     1234        ! snow heights not necessarily in height order 
     1235        !---------------------------------------------------------------- 
     1236 
     1237        DO ns = 1, m_index+1 
     1238           DO n = 0, m_index - ns + 1 
     1239              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     1240                 tmp = hitl(n) 
     1241                 hitl(n) = hitl(n+1) 
     1242                 hitl(n+1) = tmp 
     1243                 tmp = aicetl(n) 
     1244                 aicetl(n) = aicetl(n+1) 
     1245                 aicetl(n+1) = tmp 
     1246              END IF 
     1247           END DO 
     1248        END DO 
     1249 
     1250        !---------------------------------------------------------------- 
     1251        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     1252        !---------------------------------------------------------------- 
     1253 
     1254        !hitl(3)---------------------------------------------------------------- 
     1255        !                                                       | * * * * * * * * 
     1256        !                                                       |* * * * * * * * * 
     1257        !hitl(2)---------------------------------------------------------------- 
     1258        !                                    | * * * * * * * *  | * * * * * * * * 
     1259        !                                    |* * * * * * * * * |* * * * * * * * * 
     1260        !hitl(1)---------------------------------------------------------------- 
     1261        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
     1262        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
     1263        !hitl(0)---------------------------------------------------------------- 
     1264        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3) 
     1265 
     1266        ! move up over layers incrementing volume 
     1267        DO n = 1, m_index+1 
     1268 
     1269           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     1270                (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     1271 
     1272           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     1273 
     1274           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     1275              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 
     1276 
     1277              exit 
     1278           ELSE  ! still in sub-layer below the sub-layer with the depth 
     1279              rem_vol = rem_vol - vol 
     1280           END IF 
     1281 
     1282        END DO 
     1283 
     1284       END IF 
     1285 
     1286    END SUBROUTINE ice_thd_pnd_depth 
     1287 
     1288 
     1289    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 
     1290       !!------------------------------------------------------------------- 
     1291       !!                ***  ROUTINE ice_thd_pnd_perm *** 
     1292       !! 
     1293       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     1294       !!                and its permeability 
     1295       !!------------------------------------------------------------------- 
     1296 
     1297       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     1298          ticen,  &     ! internal ice temperature (K) 
     1299          salin         ! salinity (ppt)     !js: ppt according to cice 
     1300 
     1301       REAL (wp), INTENT(OUT) :: & 
     1302          perm      ! permeability 
     1303 
     1304       REAL (wp) ::   & 
     1305          Sbr       ! brine salinity 
     1306 
     1307       REAL (wp), DIMENSION(nlay_i) ::   & 
     1308          Tin, &    ! ice temperature 
     1309          phi       ! liquid fraction 
     1310 
     1311       INTEGER :: k 
     1312 
     1313       !----------------------------------------------------------------- 
     1314       ! Compute ice temperatures from enthalpies using quadratic formula 
     1315       !----------------------------------------------------------------- 
     1316 
     1317       DO k = 1,nlay_i 
     1318          Tin(k) = ticen(k) - rt0   !js: from K to degC 
     1319       END DO 
     1320 
     1321       !----------------------------------------------------------------- 
     1322       ! brine salinity and liquid fraction 
     1323       !----------------------------------------------------------------- 
     1324 
     1325       DO k = 1, nlay_i 
     1326        
     1327          ! Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 
     1328          ! Best expression to date is that one 
     1329          Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 
     1330          phi(k) = salin(k) / Sbr 
     1331           
     1332       END DO 
     1333 
     1334       !----------------------------------------------------------------- 
     1335       ! permeability 
     1336       !----------------------------------------------------------------- 
     1337 
     1338       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 
     1339 
     1340   END SUBROUTINE ice_thd_pnd_perm 
     1341    
     1342    
     1343!---------------------------------------------------------------------------------------------------------------------- 
    3041344 
    3051345   SUBROUTINE ice_thd_pnd_init  
    … …  
    3171357      INTEGER  ::   ios, ioptio   ! Local integer 
    3181358      !! 
    319       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 
    320          &                          ln_pnd_CST , rn_apnd, rn_hpnd,         & 
     1359      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 
     1360         &                          ln_pnd_CST , rn_apnd, rn_hpnd,                       & 
     1361         &                          ln_pnd_TOPO ,                                        & 
    3211362         &                          ln_pnd_lids, ln_pnd_alb 
    3221363      !!------------------------------------------------------------------- 
    … …  
    3361377         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    3371378         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd 
     1379         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO 
    3381380         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV 
    3391381         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min 
    3401382         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max 
     1383         WRITE(numout,*) '            Pond flushing efficiency                              rn_pnd_flush = ', rn_pnd_flush 
    3411384         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST 
    3421385         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd 
    … …  
    3511394      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    3521395      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF 
     1396      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF 
    3531397      IF( ioptio /= 1 )   & 
    354          & 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)' ) 
     1398         & 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)' ) 
    3551399      ! 
    3561400      SELECT CASE( nice_pnd ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_minor/src/ICE/iceupdate.F90

    r14075 r14671  
    146146            ! 
    147147            ! the non-solar is simply derived from the solar flux 
    148             qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr               
     148            qns(ji,jj) = qt_oce_ai(ji,jj) - qsr(ji,jj)               
    149149 
    150150            ! Mass flux at the atm. surface        

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

    r14075 r14671  
    244244      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
    245245      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond 
    246          &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min ) 
     246         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
    247247      END WHERE 
    248248      ! 
    … …  
    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_fix_topo_minor/src/ICE/icewri.F90

    r14075 r14671  
    104104      IF( iom_use('snwthic' ) )   CALL iom_put( 'snwthic', hm_s        * zmsk00 )                                           ! snw thickness 
    105105      IF( iom_use('icebrv'  ) )   CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 )                                           ! brine volume 
    106       IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) )          ! ice age 
     106      IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 )                                           ! ice age 
    107107      IF( iom_use('icehnew' ) )   CALL iom_put( 'icehnew', ht_i_new             )                                           ! new ice thickness formed in the leads 
    108108      IF( iom_use('snwvolu' ) )   CALL iom_put( 'snwvolu', vt_s        * zmsksn )                                           ! snow volume 
    … …  
    116116      IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth 
    117117      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area 
     118      IF( iom_use('iceepnd' ) )   CALL iom_put( 'iceepnd', SUM( a_ip_eff * a_i, dim=3 ) * zmsk00  )                         ! melt pond total effective fraction per cell area 
    118119      IF( iom_use('icehlid' ) )   CALL iom_put( 'icehlid', hm_il  * zmsk00      )                                           ! melt pond lid depth 
    119120      IF( iom_use('icevlid' ) )   CALL iom_put( 'icevlid', vt_il  * zmsk00      )                                           ! melt pond lid total volume per unit area 
    120121      ! salt 
    121       IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity 
     122      IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 )                                  ! mean ice salinity 
    122123      IF( iom_use('icesalm' ) )   CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 )                                  ! Mass of salt in sea ice per cell area 
    123124      ! heat 
    124       IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! ice mean temperature 
    125       IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) )      ! snw mean temperature 
    126       IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice surface 
    127       IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice bottom 
    128       IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the snow-ice interface 
     125      IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00  )                                      ! ice mean temperature 
     126      IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn  )                                      ! snw mean temperature 
     127      IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00  )                                      ! temperature at the ice surface 
     128      IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00  )                                      ! temperature at the ice bottom 
     129      IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00  )                                      ! temperature at the snow-ice interface 
    129130      IF( iom_use('icehc'   ) )   CALL iom_put( 'icehc'  ,  -et_i          * zmsk00 )                                       ! ice heat content 
    130131      IF( iom_use('snwhc'   ) )   CALL iom_put( 'snwhc'  ,  -et_s          * zmsksn )                                       ! snow heat content 
    … …  
    153154      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 
    154155      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories 
    155       IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories 
    156       IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories 
    157       IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories 
    158       IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age 
     156      IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l                                   ) ! thickness for categories 
     157      IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl                                   ) ! snow depth for categories 
     158      IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l                                   ) ! salinity for categories 
     159      IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l                                   ) ! ice age 
    159160      IF( iom_use('icetemp_cat' ) )   CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) & 
    160          &                                                                         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature 
     161         &                                                                         * zmsk00l                                   ) ! ice temperature 
    161162      IF( iom_use('snwtemp_cat' ) )   CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) & 
    162          &                                                                         * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature 
    163       IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature 
    164       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 
     163         &                                                                         * zmsksnl                                   ) ! snow temperature 
     164      IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l                                   ) ! surface temperature 
     165      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l                                   ) ! brine volume 
    165166      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories 
     167      IF( iom_use('icevpnd_cat' ) )   CALL iom_put( 'icevpnd_cat' ,   v_ip         * zmsk00l                                   ) ! melt pond volume for categories 
    166168      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 
    167169      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 
     170      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac per ice area for categories 
    169171      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    170172      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.