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 13809 for NEMO/branches/2020 – NEMO

Changeset 13809 for NEMO/branches/2020


Ignore:
Timestamp:
2020-11-18T11:04:09+01:00 (3 years ago)
Author:
vancop
Message:

Original melt pond code from LLN

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90

    r13284 r13809  
    3434   INTEGER ::              nice_pnd    ! choice of the type of pond scheme 
    3535   !                                   ! 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 
     36   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme 
     37   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme 
     38   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme 
     39   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme 
    3940 
    4041   !! * Substitutions 
     
    6061         ! 
    6162      CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==! 
     63         ! 
     64      CASE (np_pndTOPO)  ;   CALL pnd_TOPO   !==  Topographic melt ponds  ==! 
    6265         ! 
    6366      END SELECT 
     
    153156      !!                   Holland et al      (J. Clim, 2012) 
    154157      !!------------------------------------------------------------------- 
     158       
    155159      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    156160      !! 
     
    170174      !! 
    171175      INTEGER  ::   ji, jk                            ! loop indices 
     176       
    172177      !!------------------------------------------------------------------- 
     178       
    173179      z1_rhow   = 1._wp / rhow  
    174180      z1_aspect = 1._wp / zaspect 
     
    196202            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    197203            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?  
     204            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! >0, max for roundoff errors?  
    199205            ! 
    200206            !--- overflow ---! 
     
    208214            !    h_ip_max = 0.5 * h_i 
    209215            !    => 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) 
     216            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 
    211217            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    212218             
     
    216222            !--- Lid melting ---! 
    217223            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     224             
    218225            ! 
    219226            !--- mass flux ---! 
    220227            IF( zdv_mlt > 0._wp ) THEN 
     228               ! MV add comment on what that mass flux means 
     229               ! water removed from fw flux due to melt pond growth 
    221230               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    222231               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    223232               ! 
     233               ! MV  
     234               ! why surface melt and snow fluxes must be adjusted is not clear 
     235               ! sounds like things are counted twice 
     236               !  
    224237               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    225238               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     
    264277             
    265278            ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
     279            ! MV the expression for zsbr has wrong sign!!!! 
     280            ! MV note - the sign-corrected expression is inconsistent  
     281            ! with the rest of the SI3 code which assumes linear liquidus 
     282            ! best expression (most consistent for SI3 should be) 
     283            ! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
    266284            DO jk = 1, nlay_i 
    267285               zsbr = - 1.2_wp                                  & 
     
    269287                  &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    270288                  &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    271                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     289                
     290               ztmp(jk) = sz_i_1d(ji,jk) / zsbr ! MV -> this is brine fraction and as it reads is always negative 
    272291            END DO 
    273             zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     292            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! MV at present since ztmp is negative, this is always zero 
    274293             
    275294            ! Do the drainage using Darcy's law 
     
    301320      ! 
    302321   END SUBROUTINE pnd_LEV 
    303  
     322    
     323   SUBROUTINE pnd_TOPO    (aice,      aicen,     & 
     324                           vice,      vicen,     & 
     325                           vsnon,                & 
     326                           ticen,     salin,     & 
     327                           a_ip_frac, h_ip,      & 
     328                                      Tsfc ) 
     329                                          
     330      !!------------------------------------------------------------------- 
     331      !!                ***  ROUTINE pnd_TOPO  *** 
     332      !! 
     333      !! ** Purpose : Compute melt pond evolution 
     334      !! 
     335      !! ** Purpose :   Compute melt pond evolution based on the ice 
     336      !!                topography as inferred from the ice thickness 
     337      !!                distribution. 
     338      !! 
     339      !! ** Method  :   This code is initially based on Flocco and Feltham 
     340      !!                (2007) and Flocco et al. (2010). More to come... 
     341      !! 
     342      !! ** Tunable parameters : 
     343      !! 
     344      !! ** Note : 
     345      !! 
     346      !! ** References 
     347      !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond 
     348      !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi: 
     349      !!    10.1029/2006JC003836. 
     350      !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
     351      !!    a physically based melt pond scheme into the sea ice component of a 
     352      !!    climate model.  J. Geophys. Res. 115, C08012, 
     353      !!    doi: 10.1029/2009JC005568. 
     354      !! 
     355      !!------------------------------------------------------------------- 
     356 
     357       !js 190423: the lid on melt ponds appears only in the analog subroutine of CICE 5.1.2 
     358 
     359       REAL (wp), DIMENSION (jpi,jpj), & 
     360          INTENT(IN) :: & 
     361          aice, &    ! total ice area fraction 
     362          vice       ! total ice volume (m) 
     363 
     364       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     365          INTENT(IN) :: & 
     366          aicen, &   ! ice area fraction, per category 
     367          vsnon, &   ! snow volume, per category (m) 
     368          vicen      ! ice volume, per category (m) 
     369 
     370       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 
     371          INTENT(IN) :: & 
     372          ticen, &   ! ice temperature per category (K) 
     373          salin 
     374 
     375       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     376          INTENT(INOUT) :: & 
     377          a_ip_frac , &   ! pond area fraction of ice, per ice category 
     378          h_ip       ! pond depth, per ice category (m) 
     379 
     380       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     381          INTENT(IN) :: & 
     382          Tsfc       ! snow/sea ice surface temperature 
     383 
     384       ! local variables 
     385       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
     386          zTsfcn, & ! ice/snow surface temperature (C) 
     387          zvolpn    ! pond volume per unit area, per category (m) 
     388 
     389       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
     390          zrfrac, & ! fraction of available meltwater retained for melt ponding 
     391          zapondn,& ! pond area fraction, per category 
     392          zhpondn   ! pond depth, per category (m) 
     393 
     394       REAL (wp), DIMENSION (jpi,jpj) :: & 
     395          zvolp,    & ! total volume of pond, per unit area of pond (m) 
     396          zwfx_tmp    ! temporary array for melt water 
     397 
     398       REAL (wp) :: & 
     399          zhi,      & ! ice thickness (m) 
     400          zTavg,    & ! mean surface temperature across categories (C) 
     401          z1_rhow, & ! inverse freshwater density 
     402          zdTs,     & ! temperature difference for freeze-up (C) 
     403          zvpold,   & ! dummy pond volume 
     404          zdvn        ! change in melt pond volume for fresh water budget 
     405           
     406       INTEGER, DIMENSION (jpi*jpj) :: & 
     407          indxi, indxj    ! compressed indices for cells with ice melting 
     408 
     409       INTEGER :: ij,icells,ji,jj,jl ! loop indices 
     410 
     411       REAL (wp), PARAMETER :: & 
     412          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 
     413          zTp       = -0.15_wp, & ! pond freezing temperature (C) 
     414          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 
     415          zrexp     = 0.01_wp,  & ! constant melt pond freeze-up rate 
     416          z01       = 0.01_wp,  & 
     417          z25       = 0.25_wp,  & 
     418          z5        = 0.5_wp 
     419 
     420      z1_rhow     = 1. / rhow 
     421 
     422      !--------------------------------------------------------------- 
     423      ! Initialization 
     424      !--------------------------------------------------------------- 
     425      zhpondn  (:,:,:) = 0._wp 
     426      zapondn  (:,:,:) = 0._wp 
     427      zvolpn   (:,:,:) = 0._wp 
     428 
     429      zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0         ! Convert in Celsius 
     430 
     431      IF ( ln_pnd_fw ) THEN 
     432         v_ip_b(:,:,:) = v_ip(:,:,:) 
     433      ELSE 
     434         v_ip_b(:,:,:) = 0._wp 
     435      ENDIF 
     436 
     437      !------------------------------------------------------------------ 
     438      ! Available melt water for melt ponding and corresponding fraction 
     439      !------------------------------------------------------------------ 
     440      !js 03/05/19 unset restriction on sign of wfx_pnd_in;  mask values close to zero for future division 
     441      !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding 
     442      !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) * MAX( 0._wp, SIGN( 1._wp, wfx_sum(:,:) + wfx_snw_sum(:,:) - epsi10 ) ) 
     443      !wfx_pnd_in(:,:) = (wfx_sum(:,:) + wfx_snw_sum(:,:)) * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) 
     444 
     445      ! MV 
     446      ! NB: wfx_pnd_in can be slightly negative for very small values (why?) 
     447      ! This can in some occasions give negative 
     448      ! v_ip in the first category, which then gives crazy pond 
     449      ! fractions and crashes the code as soon as the melt-pond 
     450      ! radiative coupling is activated 
     451      ! if we understand and remove why wfx_sum or wfx_snw could be 
     452      ! negative, then, we can remove the MAX 
     453      ! NB: I now changed to wfx_snw_sum, this may fix the problem. 
     454      ! We should check 
     455 
     456 
     457      ! OLI 07/2017: when we (MV & OL) first started the inclusion of melt 
     458      ! ponds in the model, we removed the Holland et al. (2012, see 
     459      ! CESM scheme above) parameterization. I put it back here, 
     460      ! because I think it is needed. In summary, the sinks of FW for 
     461      ! ponds are: 1/ Runoff through cracks/leads => depends on the 
     462      !               total ice area only 
     463      !            2/ overflow, including Lüthje et al. (2006) limitation 
     464      !               (max a_ip fraction function of h_i). This is in 
     465      !               fact an other form of runoff that depends on the 
     466      !               ITD 
     467      !            3/ Flushing - losses by permeability 
     468      !            4/ Refreezing 
     469      !            5/ Removal of ponds on thin ice 
     470      ! I think 1 is needed because it is different from 2. However, 
     471      ! test runs should/could be done, to check the sensitivity and 
     472      ! the real usefulness of that stuff. 
     473      ! Note : the Holland et al. param was wrongly wired in NEMO3.1 (using 
     474      ! a_i instead of at_i), which might well explain why I had a too 
     475      ! weak melt pond cover in my simulations (compared to MODIS, in 
     476      ! situ obs. and CICE simulations. 
     477 
     478      !js 23/04/19: rewired back to a fraction with a_i 
     479      zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:) 
     480      zwfx_tmp(:,:) = 0._wp 
     481 
     482! MV ---> use expression from level-ice ponds, clarify what is needed 
     483 
     484      !--- Add retained melt water to melt ponds 
     485      ! v_ip should never be negative, otherwise code crashes 
     486      ! Rq MV: as far as I saw, UM5 can create very small negative v_ip values 
     487      ! hence I added the max, which was not required with Prather (1 yr run) 
     488      ! OLI: Here I use vt_ip, so I don't know if the max is 
     489      ! required... 
     490      !zvolp(:,:) = MAX(vt_ip(:,:),0._wp) + zrfrac(:,:) * wfx_pnd_in(:,:) * z1_rhow * rdt_ice  ! Total available melt water, to be distributed as melt ponds 
     491      ! OLI: 07/2017 Bugfix above, removed " * aice(:,:)" 
     492      !js 19/04/18: change zrfrac to use aicen 
     493 
     494      zvolp(:,:) = vt_ip(:,:) 
     495 
     496      DO jl = 1, jpl 
     497         ! Melt water, to be distributed as melt ponds 
     498         zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl)                                                  &  
     499                                    * ( dh_i_pnd(:,:,jl)*rhoic + dh_s_pnd(:,:,jl)*rhosn )          &  
     500                                    * z1_rhow * a_i(:,:,jl) 
     501      END DO 
     502 
     503! MV ---> use expression from level ice melt ponds (dv_mlt) 
     504 
     505      !js 03/05/19: we truncate negative values after calculating zvolp, in a 
     506      ! similar manner to the subroutine lim_mp_cesm. Variation dh_i_pnd and 
     507      ! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them 
     508      ! to be negative for some reasons. We keep this behaviour as it is, for 
     509      ! fluxes conservation reasons. If some dh are positive, then we remove water 
     510      ! indirectly from the ponds. 
     511      zvolp(:,:) = MAX( zvolp(:,:) , 0._wp ) 
     512 
     513 
     514      ! Fresh water flux going into the ponds 
     515      wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice 
     516 
     517      !--- Remove retained meltwater from surface fluxes 
     518      IF ( ln_pnd_fw ) THEN 
     519         !wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrfrac(:,:) ) 
     520         !wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrfrac(:,:) ) 
     521 
     522         !js 190419: we change the code to use a_i in zrfrac. To be tested, but 
     523         !it should be conservative. zwfx_tmp is the flux accumulated in the 
     524         !ponds. wfx_pnd_in is the total surface melt fluxes. 
     525         zwfx_tmp(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 )                           & 
     526            &            * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) 
     527         WHERE ( ABS(zwfx_tmp(:,:)) > epsi10 ) 
     528            zwfx_tmp(:,:) = wfx_pnd_in(:,:) / zwfx_tmp(:,:) 
     529         ELSEWHERE 
     530            zwfx_tmp(:,:) = 0._wp 
     531         ENDWHERE 
     532 
     533         wfx_sum(:,:)     = ( 1._wp - zwfx_tmp(:,:) ) * wfx_sum(:,:) 
     534         wfx_snw_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_snw_sum(:,:) 
     535      ENDIF 
     536 
     537       !----------------------------------------------------------------- 
     538       ! Identify grid cells with ponds 
     539       !----------------------------------------------------------------- 
     540 
     541       icells = 0 
     542       DO jj = 1, jpj 
     543         DO ji = 1, jpi 
     544            IF ( aice(ji,jj) > epsi10 ) THEN 
     545               zhi = vice(ji,jj) / aice(ji,jj) 
     546            ELSE 
     547               zhi = 0._wp 
     548            END IF 
     549 
     550            IF ( aice(ji,jj) > z01 .and. zhi > rn_himin .and. & 
     551                 zvolp(ji,jj) > zmin_volp*aice(ji,jj)) THEN 
     552               icells = icells + 1 
     553               indxi(icells) = ji 
     554               indxj(icells) = jj 
     555            ELSE  ! remove ponds on thin ice, or too small ponds 
     556               zvolpn(ji,jj,:) = 0._wp 
     557               zvolp (ji,jj) = 0._wp 
     558 
     559               a_ip(ji,jj,:)      = 0._wp 
     560               v_ip(ji,jj,:)      = 0._wp 
     561               a_ip_frac(ji,jj,:) = 0._wp 
     562               h_ip(ji,jj,:)      = 0._wp 
     563 
     564               vt_ip(ji,jj)       = 0._wp 
     565               at_ip(ji,jj)       = 0._wp 
     566 
     567!               IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean 
     568!                   wfx_pnd_out(ji,jj)   = wfx_pnd_out(ji,jj) + zvolp(ji,jj) * rhow * r1_rdtice 
     569            END IF 
     570         END DO                     ! ji 
     571      END DO                     ! jj 
     572 
     573 
     574       DO ij = 1, icells 
     575        
     576          ji = indxi(ij) 
     577          jj = indxj(ij) 
     578 
     579          !-------------------------------------------------------------- 
     580          ! Shrink pond due to refreezing 
     581          !-------------------------------------------------------------- 
     582          ! OLI 07/2017: Done like for empirical melt pond scheme (CESM). 
     583          ! Therefore, I chose to put this part of the code before the main 
     584          ! routines lim_mp_area/depth (contrary to the original code), 
     585          ! seeing the freeze-up as a global sink of 
     586          ! freshwater for melt ponds in the whole grid cell. If this was done 
     587          ! after, I would need to make an additional assumption on the shape of 
     588          ! melt ponds, which I don't want to do (for the CESM scheme, this 
     589          ! assumption was on the aspect ratio). So I remove some water due to 
     590          ! refreezing first (using zTavg instead of zTsfcn in each category) and 
     591          ! then let the FF07 routines do their job for the fractional areas and 
     592          ! depths of melt ponds. 
     593          ! The whole ice lid related stuff from FF07 was thus removed and replaced 
     594          ! by this. As mentionned below, this should be improved, but is much 
     595          ! easier to conserve heat and freshwater this way. 
     596 
     597          ! Average surface temperature is needed to compute freeze-up at the cell 
     598          ! scale 
     599          zTavg = 0._wp 
     600          DO jl = 1, jpl 
     601             zTavg = zTavg + zTsfcn(ji,jj,jl)*aicen(ji,jj,jl) 
     602          END DO 
     603          zTavg = zTavg / aice(ji,jj) 
     604 
     605          ! The freezing temperature for meltponds is assumed slightly below 0C, 
     606          ! as if meltponds had a little salt in them (hence the use of zTp). 
     607          ! The salt budget is not 
     608          ! altered for meltponds, but if it were then an actual pond freezing 
     609          ! temperature could be computed. 
     610 
     611          zdTs        = MAX ( zTp - zTavg, 0. ) 
     612 
     613          zvpold      = zvolp(ji,jj) 
     614 
     615          zvolp(ji,jj)  = zvolp(ji,jj) * EXP( zrexp * zdTs / zTp ) 
     616 
     617          !--- Dump meltwater due to refreezing ( of course this is wrong 
     618          !--- but this parameterization is too simple ) 
     619!          IF ( ln_pnd_fw ) & 
     620!             wfx_pnd_out(ji,jj)   = wfx_pnd_out(ji,jj) + rhow * ( zvpold - zvolp(ji,jj) ) * r1_rdtice 
     621!             ! OLI 07/2017 : Bugfix above, zvpold - zvolp instead of the 
     622!             ! opposite, otherwise negative contribution 
     623 
     624          !-------------------------------------------------------------- 
     625          ! calculate pond area and depth 
     626          !-------------------------------------------------------------- 
     627          zdvn = 0._wp 
     628          CALL lim_mp_area(aice(ji,jj),vice(ji,jj), & 
     629                    aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), & 
     630                    ticen(ji,jj,:,:), salin(ji,jj,:,:), & 
     631                    zvolpn(ji,jj,:), zvolp(ji,jj), & 
     632                    zapondn(ji,jj,:),zhpondn(ji,jj,:), zdvn) 
     633          ! outputs are 
     634          ! - zdvn 
     635          ! - zvolpn 
     636          ! - zvolp 
     637          ! - zapondn 
     638          ! - zhpondn 
     639 
     640!          IF ( ln_pnd_fw ) & 
     641!             wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zdvn * rhow * r1_rdtice ! update flux from ponds to ocean 
     642 
     643          !--------------------------------------------------------------- 
     644          ! Update pond volume and fraction 
     645          !--------------------------------------------------------------- 
     646          DO jl = 1, jpl 
     647             a_ip(ji,jj,jl) = zapondn(ji,jj,jl) 
     648             v_ip(ji,jj,jl) = zvolpn(ji,jj,jl) 
     649             a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / MAX(aicen(ji,jj,jl), epsi10) & 
     650                                    * MAX(0._wp, SIGN(1._wp, aicen(ji,jj,jl) - epsi10)) 
     651             h_ip     (ji,jj,jl) = zhpondn(ji,jj,jl) 
     652          END DO 
     653       END DO ! ij 
     654 
     655       IF ( ln_pnd_fw ) THEN 
     656          !js 15/05/19: water going out of the ponds give a positive freshwater 
     657          ! flux.  
     658          wfx_pnd_out(:,:) = SUM(MAX(0._wp, v_ip_b(:,:,:) - v_ip(:,:,:)), DIM=3) * rhow * r1_rdtice 
     659       ELSE 
     660         wfx_pnd_out(:,:) = 0._wp 
     661       ENDIF 
     662 
     663   END SUBROUTINE pnd_TOPO 
     664    
     665       SUBROUTINE lim_mp_area(aice,  vice, & 
     666                           aicen, vicen, vsnon, ticen, & 
     667                           salin, zvolpn, zvolp,         & 
     668                           zapondn,zhpondn,dvolp) 
     669 
     670       !!------------------------------------------------------------------- 
     671       !!                ***  ROUTINE lim_mp_area *** 
     672       !! 
     673       !! ** Purpose : Given the total volume of meltwater, update 
     674       !!              pond fraction (a_ip) and depth (should be volume) 
     675       !! 
     676       !! ** 
     677       !! 
     678       !!------------------------------------------------------------------ 
     679 
     680       REAL (wp), INTENT(IN) :: & 
     681          aice, vice 
     682 
     683       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     684          aicen, vicen, vsnon 
     685 
     686       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 
     687          ticen, salin 
     688 
     689       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 
     690          zvolpn 
     691 
     692       REAL (wp), INTENT(INOUT) :: & 
     693          zvolp, dvolp 
     694 
     695       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 
     696          zapondn, zhpondn 
     697 
     698       INTEGER :: & 
     699          n, ns,   & 
     700          m_index, & 
     701          permflag 
     702 
     703       REAL (wp), DIMENSION(jpl) :: & 
     704          hicen, & 
     705          hsnon, & 
     706          asnon, & 
     707          rhos,  &  ! OLI 07/2017 : for now this is useless, but will be useful with new snow scheme 
     708          alfan, & 
     709          betan, & 
     710          cum_max_vol, & 
     711          reduced_aicen 
     712 
     713       REAL (wp), DIMENSION(0:jpl) :: & 
     714          cum_max_vol_tmp 
     715 
     716       REAL (wp) :: & 
     717          hpond, & 
     718          drain, & 
     719          floe_weight, & 
     720          pressure_head, & 
     721          hsl_rel, & 
     722          deltah, & 
     723          perm, & 
     724          msno 
     725 
     726       REAL (wp), parameter :: & 
     727          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s 
     728          z0        = 0.0_wp    , & 
     729          c1        = 1.0_wp    , & 
     730          p4        = 0.4_wp    , & 
     731          p6        = 0.6_wp 
     732 
     733      !-----------| 
     734      !           | 
     735      !           |-----------| 
     736      !___________|___________|______________________________________sea-level 
     737      !           |           | 
     738      !           |           |---^--------| 
     739      !           |           |   |        | 
     740      !           |           |   |        |-----------|              |------- 
     741      !           |           |   |alfan(n)|           |              | 
     742      !           |           |   |        |           |--------------| 
     743      !           |           |   |        |           |              | 
     744      !---------------------------v------------------------------------------- 
     745      !           |           |   ^        |           |              | 
     746      !           |           |   |        |           |--------------| 
     747      !           |           |   |betan(n)|           |              | 
     748      !           |           |   |        |-----------|              |------- 
     749      !           |           |   |        | 
     750      !           |           |---v------- | 
     751      !           |           | 
     752      !           |-----------| 
     753      !           | 
     754      !-----------| 
     755 
     756       !------------------------------------------------------------------- 
     757       ! initialize 
     758       !------------------------------------------------------------------- 
     759 
     760       rhos(:) = rhosn ! OLI 07/2017 : same has above 
     761 
     762       DO n = 1, jpl 
     763 
     764          zapondn(n) = z0 
     765          zhpondn(n) = z0 
     766 
     767          !---------------------------------------- 
     768          ! X) compute the effective snow fraction 
     769          !---------------------------------------- 
     770          IF (aicen(n) < epsi10)  THEN 
     771             hicen(n) =  z0 
     772             hsnon(n) = z0 
     773             reduced_aicen(n) = z0 
     774             asnon(n) = z0          !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 
     775          ELSE 
     776             hicen(n) = vicen(n) / aicen(n) 
     777             hsnon(n) = vsnon(n) / aicen(n) 
     778             reduced_aicen(n) = c1 ! n=jpl 
     779 
     780             !js: initial code in NEMO_DEV 
     781             !IF (n < jpl) reduced_aicen(n) = aicen(n) & 
     782             !                     * (-0.024_wp*hicen(n) + 0.832_wp) 
     783 
     784             !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 
     785             IF (n < jpl) reduced_aicen(n) = aicen(n) & 
     786                                  * max(0.2_wp,(-0.024_wp*hicen(n) + 0.832_wp)) 
     787 
     788             asnon(n) = reduced_aicen(n)  ! effective snow fraction (empirical) 
     789             ! MV should check whether this makes sense to have the same effective snow fraction in here 
     790             ! OLI: it probably doesn't 
     791          END IF 
     792 
     793 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     794 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     795 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     796 ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     797 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     798 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
     799 ! alfan = 60% of the ice volume) in each category lies above the reference line, 
     800 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     801 
     802 ! MV: 
     803 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
     804 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
     805 
     806 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 
     807 
     808          alfan(n) = 0.6 * hicen(n) 
     809          betan(n) = 0.4 * hicen(n) 
     810 
     811          cum_max_vol(n)     = z0 
     812          cum_max_vol_tmp(n) = z0 
     813 
     814       END DO ! jpl 
     815 
     816       cum_max_vol_tmp(0) = z0 
     817       drain = z0 
     818       dvolp = z0 
     819 
     820       !---------------------------------------------------------- 
     821       ! x) Drain overflow water, update pond fraction and volume 
     822       !---------------------------------------------------------- 
     823 
     824       !-------------------------------------------------------------------------- 
     825       ! the maximum amount of water that can be contained up to each ice category 
     826       !-------------------------------------------------------------------------- 
     827 
     828       ! MV 
     829       ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     830       ! Then the excess volume cum_max_vol(jl) drains out of the system 
     831       ! It should be added to wfx_pnd_out 
     832       ! END MV 
     833       !js 18/04/19: XXX do something about this flux thing 
     834 
     835       DO n = 1, jpl-1 ! last category can not hold any volume 
     836 
     837          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 
     838 
     839             ! total volume in level including snow 
     840             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 
     841                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 
     842 
     843             ! subtract snow solid volumes from lower categories in current level 
     844             DO ns = 1, n 
     845                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 
     846                   - rhos(ns)/rhow * &   ! free air fraction that can be filled by water 
     847                     asnon(ns)  * &    ! effective areal fraction of snow in that category 
     848                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), z0) 
     849             END DO 
     850 
     851          ELSE ! assume higher categories unoccupied 
     852             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 
     853          END IF 
     854          !IF (cum_max_vol_tmp(n) < z0) THEN 
     855          !   CALL abort_ice('negative melt pond volume') 
     856          !END IF 
     857       END DO 
     858       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     859       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     860 
     861       !---------------------------------------------------------------- 
     862       ! is there more meltwater than can be held in the floe? 
     863       !---------------------------------------------------------------- 
     864       IF (zvolp >= cum_max_vol(jpl)) THEN 
     865          drain = zvolp - cum_max_vol(jpl) + epsi10 
     866          zvolp = zvolp - drain ! update meltwater volume available 
     867          dvolp = drain         ! this is the drained water 
     868          IF (zvolp < epsi10) THEN 
     869             dvolp = dvolp + zvolp 
     870             zvolp = z0 
     871          END IF 
     872       END IF 
     873 
     874       ! height and area corresponding to the remaining volume 
     875 
     876      CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 
     877 
     878       DO n=1, m_index 
     879          !zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 
     880          !                                         !  volume instead, no ? 
     881          zhpondn(n) = max((hpond - alfan(n) + alfan(1)), z0)      !js: from CICE 5.1.2 
     882          zapondn(n) = reduced_aicen(n) 
     883          ! in practise, pond fraction depends on the empirical snow fraction 
     884          ! so in turn on ice thickness 
     885       END DO 
     886       !zapond = sum(zapondn(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     887 
     888       !------------------------------------------------------------------------ 
     889       ! Drainage through brine network (permeability) 
     890       !------------------------------------------------------------------------ 
     891       !!! drainage due to ice permeability - Darcy's law 
     892 
     893       ! sea water level 
     894       msno = z0 
     895       DO n=1,jpl 
     896         msno = msno + vsnon(n) * rhos(n) 
     897       END DO 
     898       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 
     899       hsl_rel = floe_weight / rau0 & 
     900               - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 
     901 
     902       deltah = hpond - hsl_rel 
     903       pressure_head = grav * rau0 * max(deltah, z0) 
     904 
     905       ! drain if ice is permeable 
     906       permflag = 0 
     907       IF (pressure_head > z0) THEN 
     908          DO n = 1, jpl-1 
     909             IF (hicen(n) /= z0) THEN 
     910             !IF (hicen(n) > z0) THEN           !js: from CICE 5.1.2 
     911                perm = 0._wp ! MV ugly dummy patch 
     912                CALL lim_mp_perm(ticen(:,n), salin(:,n), perm) 
     913                IF (perm > z0) permflag = 1 
     914 
     915                drain = perm*zapondn(n)*pressure_head*rdt_ice / & 
     916                                         (viscosity*hicen(n)) 
     917                dvolp = dvolp + min(drain, zvolp) 
     918                zvolp = max(zvolp - drain, z0) 
     919                IF (zvolp < epsi10) THEN 
     920                   dvolp = dvolp + zvolp 
     921                   zvolp = z0 
     922                END IF 
     923            END IF 
     924         END DO 
     925 
     926          ! adjust melt pond dimensions 
     927          IF (permflag > 0) THEN 
     928             ! recompute pond depth 
     929            CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 
     930             DO n=1, m_index 
     931                zhpondn(n) = hpond - alfan(n) + alfan(1) 
     932                zapondn(n) = reduced_aicen(n) 
     933             END DO 
     934             !zapond = sum(zapondn(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     935          END IF 
     936       END IF ! pressure_head 
     937 
     938       !------------------------------- 
     939       ! X) remove water from the snow 
     940       !------------------------------- 
     941       !------------------------------------------------------------------------ 
     942       ! total melt pond volume in category does not include snow volume 
     943       ! snow in melt ponds is not melted 
     944       !------------------------------------------------------------------------ 
     945 
     946       ! Calculate pond volume for lower categories 
     947       DO n=1,m_index-1 
     948          zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 
     949                   - (rhos(n)/rhow) * asnon(n) * min(hsnon(n), zhpondn(n)) 
     950       END DO 
     951 
     952       ! Calculate pond volume for highest category = remaining pond volume 
     953 
     954       ! The following is completely unclear to Martin at least 
     955       ! Could we redefine properly and recode in a more readable way ? 
     956 
     957       ! m_index = last category with melt pond 
     958 
     959       IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 
     960 
     961       IF (m_index > 1) THEN 
     962         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 
     963           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) 
     964         ELSE 
     965           zvolpn(m_index) = z0 
     966           zhpondn(m_index) = z0 
     967           zapondn(m_index) = z0 
     968           ! If remaining pond volume is negative reduce pond volume of 
     969           ! lower category 
     970           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 
     971             zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp 
     972         END IF 
     973       END IF 
     974 
     975       DO n=1,m_index 
     976          IF (zapondn(n) > epsi10) THEN 
     977              zhpondn(n) = zvolpn(n) / zapondn(n) 
     978          ELSE 
     979             dvolp = dvolp + zvolpn(n) 
     980             zhpondn(n) = z0 
     981             zvolpn(n) = z0 
     982             zapondn(n) = z0 
     983          END IF 
     984       END DO 
     985       DO n = m_index+1, jpl 
     986          zhpondn(n) = z0 
     987          zapondn(n) = z0 
     988          zvolpn (n) = z0 
     989       END DO 
     990 
     991    END SUBROUTINE lim_mp_area 
     992 
     993 
     994    SUBROUTINE lim_mp_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 
     995       !!------------------------------------------------------------------- 
     996       !!                ***  ROUTINE lim_mp_depth  *** 
     997       !! 
     998       !! ** Purpose :   Compute melt pond depth 
     999       !!------------------------------------------------------------------- 
     1000 
     1001       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     1002          aicen, & 
     1003          asnon, & 
     1004          hsnon, & 
     1005          rhos,  & 
     1006          alfan, & 
     1007          cum_max_vol 
     1008 
     1009       REAL (wp), INTENT(IN) :: & 
     1010          zvolp 
     1011 
     1012       REAL (wp), INTENT(OUT) :: & 
     1013          hpond 
     1014 
     1015       INTEGER, INTENT(OUT) :: & 
     1016          m_index 
     1017 
     1018       INTEGER :: n, ns 
     1019 
     1020       REAL (wp), DIMENSION(0:jpl+1) :: & 
     1021          hitl, & 
     1022          aicetl 
     1023 
     1024       REAL (wp) :: & 
     1025          rem_vol, & 
     1026          area, & 
     1027          vol, & 
     1028          tmp, & 
     1029          z0   = 0.0_wp 
     1030 
     1031       !---------------------------------------------------------------- 
     1032       ! hpond is zero if zvolp is zero - have we fully drained? 
     1033       !---------------------------------------------------------------- 
     1034 
     1035       IF (zvolp < epsi10) THEN 
     1036        hpond = z0 
     1037        m_index = 0 
     1038       ELSE 
     1039 
     1040        !---------------------------------------------------------------- 
     1041        ! Calculate the category where water fills up to 
     1042        !---------------------------------------------------------------- 
     1043 
     1044        !----------| 
     1045        !          | 
     1046        !          | 
     1047        !          |----------|                                     -- -- 
     1048        !__________|__________|_________________________________________ ^ 
     1049        !          |          |             rem_vol     ^                | Semi-filled 
     1050        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     1051        !          |          |          |              | 
     1052        !          |          |          |              |hpond 
     1053        !          |          |          |----------|   |     |------- 
     1054        !          |          |          |          |   |     | 
     1055        !          |          |          |          |---v-----| 
     1056        !          |          | m_index  |          |         | 
     1057        !------------------------------------------------------------- 
     1058 
     1059        m_index = 0  ! 1:m_index categories have water in them 
     1060        DO n = 1, jpl 
     1061           IF (zvolp <= cum_max_vol(n)) THEN 
     1062              m_index = n 
     1063              IF (n == 1) THEN 
     1064                 rem_vol = zvolp 
     1065              ELSE 
     1066                 rem_vol = zvolp - cum_max_vol(n-1) 
     1067              END IF 
     1068              exit ! to break out of the loop 
     1069           END IF 
     1070        END DO 
     1071        m_index = min(jpl-1, m_index) 
     1072 
     1073        !---------------------------------------------------------------- 
     1074        ! semi-filled layer may have m_index different snow in it 
     1075        !---------------------------------------------------------------- 
     1076 
     1077        !-----------------------------------------------------------  ^ 
     1078        !                                                             |  alfan(m_index+1) 
     1079        !                                                             | 
     1080        !hitl(3)-->                             |----------|          | 
     1081        !hitl(2)-->                |------------| * * * * *|          | 
     1082        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     1083        !hitl(0)-->-------------------------------------------------  |  ^ 
     1084        !                various snow from lower categories          |  |alfa(m_index) 
     1085 
     1086        ! hitl - heights of the snow layers from thinner and current categories 
     1087        ! aicetl - area of each snow depth in this layer 
     1088 
     1089        hitl(:) = z0 
     1090        aicetl(:) = z0 
     1091        DO n = 1, m_index 
     1092           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     1093                                  alfan(m_index+1) - alfan(m_index)), z0) 
     1094           aicetl(n) = asnon(n) 
     1095 
     1096           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     1097        END DO 
     1098 
     1099        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     1100        aicetl(m_index+1) = z0 
     1101 
     1102        !---------------------------------------------------------------- 
     1103        ! reorder array according to hitl 
     1104        ! snow heights not necessarily in height order 
     1105        !---------------------------------------------------------------- 
     1106 
     1107        DO ns = 1, m_index+1 
     1108           DO n = 0, m_index - ns + 1 
     1109              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     1110                 tmp = hitl(n) 
     1111                 hitl(n) = hitl(n+1) 
     1112                 hitl(n+1) = tmp 
     1113                 tmp = aicetl(n) 
     1114                 aicetl(n) = aicetl(n+1) 
     1115                 aicetl(n+1) = tmp 
     1116              END IF 
     1117           END DO 
     1118        END DO 
     1119 
     1120        !---------------------------------------------------------------- 
     1121        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     1122        !---------------------------------------------------------------- 
     1123 
     1124        !hitl(3)---------------------------------------------------------------- 
     1125        !                                                       | * * * * * * * * 
     1126        !                                                       |* * * * * * * * * 
     1127        !hitl(2)---------------------------------------------------------------- 
     1128        !                                    | * * * * * * * *  | * * * * * * * * 
     1129        !                                    |* * * * * * * * * |* * * * * * * * * 
     1130        !hitl(1)---------------------------------------------------------------- 
     1131        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
     1132        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
     1133        !hitl(0)---------------------------------------------------------------- 
     1134        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3) 
     1135 
     1136        ! move up over layers incrementing volume 
     1137        DO n = 1, m_index+1 
     1138 
     1139           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     1140                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     1141 
     1142           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     1143 
     1144           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     1145              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 
     1146 
     1147              exit 
     1148           ELSE  ! still in sub-layer below the sub-layer with the depth 
     1149              rem_vol = rem_vol - vol 
     1150           END IF 
     1151 
     1152        END DO 
     1153 
     1154       END IF 
     1155 
     1156    END SUBROUTINE lim_mp_depth 
     1157 
     1158 
     1159    SUBROUTINE lim_mp_perm(ticen, salin, perm) 
     1160       !!------------------------------------------------------------------- 
     1161       !!                ***  ROUTINE lim_mp_perm *** 
     1162       !! 
     1163       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     1164       !!                and its permeability 
     1165       !!------------------------------------------------------------------- 
     1166 
     1167       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     1168          ticen,  &     ! internal ice temperature (K) 
     1169          salin         ! salinity (ppt)     !js: ppt according to cice 
     1170 
     1171       REAL (wp), INTENT(OUT) :: & 
     1172          perm      ! permeability 
     1173 
     1174       REAL (wp) ::   & 
     1175          Sbr       ! brine salinity 
     1176 
     1177       REAL (wp), DIMENSION(nlay_i) ::   & 
     1178          Tin, &    ! ice temperature 
     1179          phi       ! liquid fraction 
     1180 
     1181       INTEGER :: k 
     1182 
     1183       !----------------------------------------------------------------- 
     1184       ! Compute ice temperatures from enthalpies using quadratic formula 
     1185       !----------------------------------------------------------------- 
     1186 
     1187       DO k = 1,nlay_i 
     1188          Tin(k) = ticen(k) - rt0   !js: from K to degC 
     1189       END DO 
     1190 
     1191       !----------------------------------------------------------------- 
     1192       ! brine salinity and liquid fraction 
     1193       !----------------------------------------------------------------- 
     1194 
     1195       IF (maxval(Tin) <= -2.0_wp) THEN 
     1196 
     1197          ! Assur 1958 
     1198          DO k = 1,nlay_i 
     1199             Sbr = - 1.2_wp                 & 
     1200                   - 21.8_wp    * Tin(k)    & 
     1201                   - 0.919_wp   * Tin(k)**2 & 
     1202                   - 0.01878_wp * Tin(k)**3 
     1203             phi(k) = salin(k) / Sbr ! liquid fraction 
     1204 
     1205             !js: Sbr must not be zero 
     1206         END DO ! k 
     1207 
     1208       ELSE 
     1209 
     1210          ! Notz 2005 thesis eq. 3.2 
     1211          !js: "interstitial brine Sbr (in ppt) as a function of temperature T (in degC)". 
     1212          DO k = 1,nlay_i 
     1213             Sbr = - 17.6_wp   * Tin(k)    & 
     1214                   - 0.389_wp  * Tin(k)**2 & 
     1215                   - 0.00362_wp* Tin(k)**3 
     1216             phi(k) = salin(k) / Sbr ! liquid fraction 
     1217 
     1218             !js: Sbr must not be zero 
     1219          END DO 
     1220 
     1221       END IF 
     1222 
     1223       !----------------------------------------------------------------- 
     1224       ! permeability 
     1225       !----------------------------------------------------------------- 
     1226 
     1227       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 
     1228 
     1229   END SUBROUTINE lim_mp_perm 
     1230    
     1231    
     1232!---------------------------------------------------------------------------------------------------------------------- 
    3041233 
    3051234   SUBROUTINE ice_thd_pnd_init  
     
    3191248      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 
    3201249         &                          ln_pnd_CST , rn_apnd, rn_hpnd,         & 
     1250         &                          ln_pnd_TOPO ,                          & 
    3211251         &                          ln_pnd_lids, ln_pnd_alb 
    3221252      !!------------------------------------------------------------------- 
     
    3361266         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    3371267         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd 
     1268         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO 
    3381269         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV 
    3391270         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min 
Note: See TracChangeset for help on using the changeset viewer.