Changeset 7325


Ignore:
Timestamp:
2016-11-23T15:50:17+01:00 (5 years ago)
Author:
vancop
Message:

Unfinished work toward incorporation of melt ponds

Location:
branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90

    r7293 r7325  
    1414   !!   lim_mp_init      : some initialization and namelist read 
    1515   !!   lim_mp           : main calling routine 
    16    !!   compute_mp_topo  : actual melt pond routine 
    17    !!   pond_area        : computes melt pond fraction per category 
     16   !!   lim_mp_topo      : main melt pond routine for the "topographic" formulation (FloccoFeltham) 
     17   !!   lim_mp_area      : ??? compute melt pond fraction per category 
     18   !!   lim_mp_perm      : computes permeability (should be a FUNCTION!) 
    1819   !!   calc_hpond       : computes melt pond depth 
    1920   !!   permeability_phy : computes permeability 
    2021 
    2122   !!---------------------------------------------------------------------- 
    22 !  USE phycst           ! physical constants 
    23 !  USE dom_oce          ! ocean space and time domain 
     23   USE phycst           ! physical constants 
     24   USE dom_oce          ! ocean space and time domain 
    2425!  USE sbc_ice          ! Surface boundary condition: ice   fields 
    2526   USE ice              ! LIM-3 variables 
     
    5859 
    5960   PUBLIC   lim_mp_init    ! routine called by sbcice_lim.F90 
    60 !  PUBLIC   lim_mp         ! routine called by sbcice_lim.F90 
     61   PUBLIC   lim_mp         ! routine called by sbcice_lim.F90 
    6162 
    6263   !! * Substitutions 
     
    104105   END SUBROUTINE lim_mp_init 
    105106 
    106 #else 
    107    !!---------------------------------------------------------------------- 
    108    !!   Default option          Empty module           NO LIM sea-ice model 
    109    !!---------------------------------------------------------------------- 
    110 CONTAINS 
    111    SUBROUTINE lim_mp_init     ! Empty routine 
    112    END SUBROUTINE lim_mp_init 
    113 #endif  
    114  
    115    !!====================================================================== 
    116 END MODULE limmp  
    117  
    118 !====================================================================================================================================================== 
    119 !====================================================================================================================================================== 
    120 !====================================================================================================================================================== 
    121 ! DIRTY LEFT OVERS FROM ASSHOLES 
    122 ! 
    123 !  SUBROUTINE lim_dyn( kt ) 
    124 !     !!------------------------------------------------------------------- 
    125 !     !!               ***  ROUTINE lim_dyn  *** 
    126 !     !!                
    127 !     !! ** Purpose :   compute ice velocity 
    128 !     !!                 
    129 !     !! ** Method  :  
    130 !     !! 
    131 !     !! ** Action  : - Initialisation 
    132 !     !!              - Call of the dynamic routine for each hemisphere 
    133 !     !!------------------------------------------------------------------------------------ 
    134 !     INTEGER, INTENT(in) ::   kt     ! number of iteration 
    135 !     !! 
    136 !     INTEGER  :: jl, jk ! dummy loop indices 
    137 !     REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    138 !    !!--------------------------------------------------------------------- 
    139  
    140 !     IF( nn_timing == 1 )  CALL timing_start('limdyn') 
    141  
    142 !     CALL lim_var_agg(1)                      ! aggregate ice categories 
     107   SUBROUTINE lim_mp( kt ) 
     108      !!------------------------------------------------------------------- 
     109      !!               ***  ROUTINE lim_mp   *** 
     110      !!                
     111      !! ** Purpose :   change melt pond fraction 
     112      !!                 
     113      !! ** Method  :   brutal force 
     114      !! 
     115      !! ** Action  : -  
     116      !!              -  
     117      !!------------------------------------------------------------------------------------ 
     118 
     119      INTEGER, INTENT(in) :: kt    ! number of iteration 
     120      INTEGER  ::   ji, jj, jl     ! dummy loop indices 
     121 
     122!     REAL(wp), POINTER, DIMENSION(:,:) ::   zfsurf  ! surface heat flux(obsolete, should be droped) 
    143123      ! 
    144 !     ! conservation test 
    145 !     IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    146 !  END SUBROUTINE lim_dyn 
     124      !!------------------------------------------------------------------- 
     125 
     126      IF( nn_timing == 1 )  CALL timing_start('limthd') 
     127 
     128      IF( nn_timing == 1 )  CALL timing_start('lim_mp') 
     129 
     130      CALL lim_mp_topo    (at_i, a_i,                                       & 
     131                &          vt_i, v_i, v_s,            t_i, s_i, a_ip_frac,  & 
     132                &          h_ip,     t_su) 
     133 
     134      ! we should probably not aggregate here since we do it in lim_var_agg 
     135      ! before output, unless we need the total volume and faction else where 
     136 
     137      ! we should also make sure a_ip and v_ip are properly updated at the end 
     138      ! of the routine 
     139 
     140   END SUBROUTINE lim_mp  
     141 
     142   SUBROUTINE lim_mp_topo    (aice,      aicen,     & 
     143                              vice,      vicen,     & 
     144                              vsnon,                & 
     145                              ticen,     salin,     & 
     146                              a_ip_frac, h_ip,      & 
     147                                         Tsfc ) 
     148       !!------------------------------------------------------------------- 
     149       !!                ***  ROUTINE lim_mp_topo  *** 
     150       !! 
     151       !! ** Purpose :   Compute melt pond evolution based on the ice  
     152       !!                topography as inferred from the ice thickness  
     153       !!                distribution.   
     154       !! 
     155       !! ** Method  :   This code is initially based on Flocco and Feltham  
     156       !!                (2007) and Flocco et al. (2010). More to come... 
     157       !! 
     158       !! ** Tunable parameters : 
     159       !!  
     160       !! ** Note : 
     161       !! 
     162       !! ** References 
     163       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond  
     164       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:  
     165       !!    10.1029/2006JC003836. 
     166       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
     167       !!    a physically based melt pond scheme into the sea ice component of a 
     168       !!    climate model.  J. Geophys. Res. 115, C08012,  
     169       !!    doi: 10.1029/2009JC005568. 
     170       !!     
     171       !!------------------------------------------------------------------- 
     172  
     173       REAL (wp), DIMENSION (jpi,jpj), & 
     174          INTENT(IN) :: & 
     175          aice, &    ! total ice area fraction 
     176          vice       ! total ice volume (m) 
     177  
     178       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     179          INTENT(IN) :: & 
     180          aicen, &   ! ice area fraction, per category 
     181          vsnon, &   ! snow volume, per category (m) 
     182          vicen      ! ice volume, per category (m) 
     183  
     184       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 
     185          INTENT(IN) :: & 
     186          ticen, &   ! ice enthalpy, per category 
     187          salin 
     188  
     189       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     190          INTENT(INOUT) :: & 
     191          a_ip_frac , &   ! pond area fraction of ice, per ice category 
     192          h_ip       ! pond depth, per ice category (m) 
     193  
     194       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
     195          INTENT(IN) :: & 
     196          Tsfc       ! snow/sea ice surface temperature 
     197  
     198       ! local variables 
     199       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
     200          zTsfcn, & ! ice/snow surface temperature (C) 
     201          zvolpn, & ! pond volume per unit area, per category (m) 
     202          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m) 
     203  
     204       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
     205          zapondn,& ! pond area fraction, per category 
     206          zhpondn   ! pond depth, per category (m) 
     207  
     208       REAL (wp), DIMENSION (jpi,jpj) :: & 
     209          zvolp       ! total volume of pond, per unit area of pond (m) 
     210  
     211       REAL (wp) :: & 
     212          zhi,    & ! ice thickness (m) 
     213          zdHui,  & ! change in thickness of ice lid (m) 
     214          zomega, & ! conduction 
     215          zdTice, & ! temperature difference across ice lid (C) 
     216          zdvice, & ! change in ice volume (m) 
     217          zTavg,  & ! mean surface temperature across categories (C) 
     218          zTp,    & ! pond freezing temperature (C) 
     219          zdvn      ! change in melt pond volume for fresh water budget 
     220       INTEGER, DIMENSION (jpi*jpj) :: & 
     221          indxi, indxj    ! compressed indices for cells with ice melting 
     222  
     223       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices 
     224  
     225       INTEGER, DIMENSION (jpl) :: & 
     226          kcells          ! cells where ice lid combines with vice 
     227  
     228       INTEGER, DIMENSION (jpi*jpj,jpl) :: & 
     229          indxii, indxjj  ! i,j indices for kcells loop 
     230  
     231       REAL (wp), parameter :: & 
     232          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m)  
     233          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C) 
     234          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 
     235          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 
     236          z0       = 0._wp,    &  
     237          zTimelt   = 0._wp,    & 
     238          z01      = 0.01_wp,  & 
     239          z25      = 0.25_wp,  & 
     240          z5       = 0.5_wp 
     241 
     242       !--------------------------------------------------------------- 
     243       ! Initialization 
     244       !--------------------------------------------------------------- 
     245  
     246       zhpondn(:,:,:) = 0._wp 
     247       zapondn(:,:,:) = 0._wp 
     248       indxii(:,:) = 0 
     249       indxjj(:,:) = 0 
     250       kcells(:)   = 0 
     251 
     252       zvolp(:,:) = wfx_sum(:,:) + wfx_snw(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds 
     253       zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0                   ! Convert in Celsius 
     254  
     255       ! The freezing temperature for meltponds is assumed slightly below 0C, 
     256       ! as if meltponds had a little salt in them.  The salt budget is not 
     257       ! altered for meltponds, but if it were then an actual pond freezing  
     258       ! temperature could be computed. 
     259  
     260       ! zTp = zTimelt - zTd  ---> for lids 
     261  
     262       !----------------------------------------------------------------- 
     263       ! Identify grid cells with ponds 
     264       !----------------------------------------------------------------- 
     265  
     266       icells = 0 
     267       DO j = 1, jpj 
     268       DO i = 1, jpi 
     269          zhi = z0 
     270          IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j) 
     271          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 
     272             zvolp(i,j) > zmin_volp*aice(i,j)) THEN 
     273             icells = icells + 1 
     274             indxi(icells) = i 
     275             indxj(icells) = j 
     276          ELSE  ! remove ponds on thin ice 
     277             !fpond(i,j) = fpond(i,j) - zvolp(i,j) 
     278             zvolpn(i,j,:) = z0 
     279             zvuin (i,j,:) = z0 
     280             zvolp (i,j) = z0 
     281          END IF 
     282       END DO                     ! i 
     283       END DO                     ! j 
     284  
     285       DO ij = 1, icells 
     286          i = indxi(ij) 
     287          j = indxj(ij) 
     288  
     289          !-------------------------------------------------------------- 
     290          ! calculate pond area and depth 
     291          !-------------------------------------------------------------- 
     292          CALL lim_mp_area(aice(i,j),vice(i,j), & 
     293                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), & 
     294                    ticen(i,j,:,:), salin(i,j,:,:), & 
     295                    zvolpn(i,j,:), zvolp(i,j), & 
     296                    zapondn(i,j,:),zhpondn(i,j,:), zdvn) 
     297         ! outputs are 
     298         ! - zdvn 
     299         ! - zvolpn 
     300         ! - zvolp 
     301         ! - zapondn 
     302         ! - zhpondn 
     303 
     304          wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean 
     305  
     306          ! mean surface temperature 
     307          zTavg = z0 
     308          DO n = 1, jpl 
     309             zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n) 
     310          END DO 
     311          zTavg = zTavg / aice(i,j) 
     312  
     313       END DO ! ij 
     314  
     315       !--------------------------------------------------------------- 
     316       ! Update pond volume and fraction 
     317       !--------------------------------------------------------------- 
     318  
     319       a_ip(:,:,:) = zapondn(:,:,:) 
     320       v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:) 
     321       a_ip_frac(:,:,:) = 0._wp 
     322       h_ip     (:,:,:) = 0._wp 
     323 
     324    END SUBROUTINE lim_mp_topo 
     325 
     326    SUBROUTINE lim_mp_area(aice,vice, & 
     327                         aicen, vicen, vsnon, ticen, & 
     328                         salin, zvolpn, zvolp,         & 
     329                         zapondn,zhpondn,dvolp) 
     330 
     331       !!------------------------------------------------------------------- 
     332       !!                ***  ROUTINE lim_mp_area *** 
     333       !! 
     334       !! ** Purpose : Given the total volume of meltwater, update 
     335       !!              pond fraction (a_ip) and depth (should be volume) 
     336       !! 
     337       !! ** 
     338       !!               
     339       !!------------------------------------------------------------------ 
     340 
     341       REAL (wp), INTENT(IN) :: & 
     342          aice,vice 
     343  
     344       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     345          aicen, vicen, vsnon 
     346  
     347       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 
     348          ticen, salin 
     349  
     350       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 
     351          zvolpn 
     352  
     353       REAL (wp), INTENT(INOUT) :: & 
     354          zvolp, dvolp 
     355  
     356       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 
     357          zapondn, zhpondn 
     358  
     359       INTEGER :: & 
     360          n, ns,   & 
     361          m_index, & 
     362          permflag 
     363  
     364       REAL (wp), DIMENSION(jpl) :: & 
     365          hicen, & 
     366          hsnon, & 
     367          asnon, & 
     368          alfan, & 
     369          betan, & 
     370          cum_max_vol, & 
     371          reduced_aicen         
     372  
     373       REAL (wp), DIMENSION(0:jpl) :: & 
     374          cum_max_vol_tmp 
     375  
     376       REAL (wp) :: & 
     377          hpond, & 
     378          drain, & 
     379          floe_weight, & 
     380          pressure_head, & 
     381          hsl_rel, & 
     382          deltah, & 
     383          perm, & 
     384          msno 
     385  
     386       REAL (wp), parameter :: &  
     387          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s 
     388          z0        = 0.0_wp    , & 
     389          c1        = 1.0_wp    , & 
     390          p4        = 0.4_wp    , & 
     391          p6        = 0.6_wp    , & 
     392          epsi10      = 1.0e-11_wp 
     393           
     394      !-----------| 
     395      !           | 
     396      !           |-----------| 
     397      !___________|___________|______________________________________sea-level 
     398      !           |           | 
     399      !           |           |---^--------| 
     400      !           |           |   |        | 
     401      !           |           |   |        |-----------|              |------- 
     402      !           |           |   |alfan(n)|           |              | 
     403      !           |           |   |        |           |--------------| 
     404      !           |           |   |        |           |              | 
     405      !---------------------------v------------------------------------------- 
     406      !           |           |   ^        |           |              | 
     407      !           |           |   |        |           |--------------| 
     408      !           |           |   |betan(n)|           |              | 
     409      !           |           |   |        |-----------|              |------- 
     410      !           |           |   |        | 
     411      !           |           |---v------- | 
     412      !           |           | 
     413      !           |-----------| 
     414      !           | 
     415      !-----------| 
     416      
     417       !------------------------------------------------------------------- 
     418       ! initialize 
     419       !------------------------------------------------------------------- 
     420  
     421       DO n = 1, jpl 
     422  
     423          zapondn(n) = z0 
     424          zhpondn(n) = z0 
     425  
     426          !---------------------------------------- 
     427          ! X) compute the effective snow fraction 
     428          !---------------------------------------- 
     429          IF (aicen(n) < epsi10)  THEN 
     430             hicen(n) =  z0  
     431             hsnon(n) = z0 
     432             reduced_aicen(n) = z0 
     433          ELSE 
     434             hicen(n) = vicen(n) / aicen(n) 
     435             hsnon(n) = vsnon(n) / aicen(n) 
     436             reduced_aicen(n) = c1 ! n=jpl 
     437             IF (n < jpl) reduced_aicen(n) = aicen(n) & 
     438                                  * (-0.024_wp*hicen(n) + 0.832_wp)  
     439             asnon(n) = reduced_aicen(n)  ! effective snow fraction (empirical) 
     440          END IF 
     441  
     442 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     443 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     444 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     445 ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     446 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     447 ! Here, alfa = 60% of the ice area (and since hice is constant in a category,  
     448 ! alfan = 60% of the ice volume) in each category lies above the reference line,  
     449 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     450  
     451          alfan(n) = p6 * hicen(n) 
     452          betan(n) = p4 * hicen(n) 
     453         
     454          cum_max_vol(n)     = z0 
     455          cum_max_vol_tmp(n) = z0 
     456      
     457       END DO ! jpl 
     458  
     459       cum_max_vol_tmp(0) = z0 
     460       drain = z0 
     461       dvolp = z0 
     462 
     463       !---------------------------------------------------------- 
     464       ! x) Drain overflow water, update pond fraction and volume 
     465       !---------------------------------------------------------- 
     466        
     467       !-------------------------------------------------------------------------- 
     468       ! the maximum amount of water that can be contained up to each ice category 
     469       !-------------------------------------------------------------------------- 
     470 
     471       ! MV 
     472       ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     473       ! Then the excess volume cum_max_vol(jl) drains out of the system 
     474       ! It should be added to wfx_pnd 
     475       ! END MV 
     476      
     477       DO n = 1, jpl-1 ! last category can not hold any volume 
     478  
     479          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 
     480  
     481             ! total volume in level including snow 
     482             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 
     483                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))  
     484  
     485             ! subtract snow solid volumes from lower categories in current level 
     486             DO ns = 1, n  
     487                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 
     488                   - rhosn/rhofw * &   ! free air fraction that can be filled by water 
     489                     asnon(ns)  * &    ! effective areal fraction of snow in that category 
     490                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- & 
     491                                   alfan(n)), z0)   
     492             END DO 
     493 
     494          ELSE ! assume higher categories unoccupied 
     495             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 
     496          END IF 
     497          !IF (cum_max_vol_tmp(n) < z0) THEN 
     498          !   call abort_ice('negative melt pond volume') 
     499          !END IF 
     500       END DO 
     501       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     502       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     503      
     504       !---------------------------------------------------------------- 
     505       ! is there more meltwater than can be held in the floe? 
     506       !---------------------------------------------------------------- 
     507       IF (zvolp >= cum_max_vol(jpl)) THEN 
     508          drain = zvolp - cum_max_vol(jpl) + epsi10 
     509          zvolp = zvolp - drain ! update meltwater volume available 
     510          dvolp = drain         ! this is the drained water 
     511          IF (zvolp < epsi10) THEN 
     512             dvolp = dvolp + zvolp 
     513             zvolp = z0 
     514          END IF 
     515       END IF 
     516      
     517       ! height and area corresponding to the remaining volume 
     518  
     519!      call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
     520!           zvolp, cum_max_vol, hpond, m_index) 
     521      
     522       DO n=1, m_index 
     523          zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 
     524                                                   !  volume instead, no ? 
     525          zapondn(n) = reduced_aicen(n)  
     526          ! in practise, pond fraction depends on the empirical snow fraction 
     527          ! so in turn on ice thickness 
     528       END DO 
     529      
     530       !------------------------------------------------------------------------ 
     531       ! Drainage through brine network (permeability) 
     532       !------------------------------------------------------------------------ 
     533       !!! drainage due to ice permeability - Darcy's law 
     534      
     535       ! sea water level  
     536       msno = z0 
     537       DO n=1,jpl 
     538         msno = msno + vsnon(n) * rhosn 
     539       END DO 
     540       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 
     541       hsl_rel = floe_weight / rau0 & 
     542               - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 
     543      
     544       deltah = hpond - hsl_rel 
     545       pressure_head = grav * rau0 * max(deltah, z0) 
     546  
     547       ! drain IF ice is permeable     
     548       permflag = 0 
     549       IF (pressure_head > z0) THEN 
     550       DO n = 1, jpl-1 
     551          IF (hicen(n) /= z0) THEN 
     552             perm = 0. ! MV ugly dummy patch 
     553             CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm) 
     554             IF (perm > z0) permflag = 1 
     555 
     556             drain = perm*zapondn(n)*pressure_head*rdt_ice / & 
     557                                      (viscosity*hicen(n)) 
     558             dvolp = dvolp + min(drain, zvolp) 
     559             zvolp = max(zvolp - drain, z0) 
     560             IF (zvolp < epsi10) THEN 
     561                dvolp = dvolp + zvolp 
     562                zvolp = z0 
     563             END IF 
     564          END IF 
     565       END DO 
     566   
     567       ! adjust melt pond DIMENSIONs 
     568       IF (permflag > 0) THEN 
     569          ! recompute pond depth     
     570!         CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
     571!                         zvolp, cum_max_vol, hpond, m_index) 
     572          DO n=1, m_index 
     573             zhpondn(n) = hpond - alfan(n) + alfan(1) 
     574             zapondn(n) = reduced_aicen(n)  
     575          END DO 
     576       END IF 
     577       END IF ! pressure_head 
     578  
     579       !------------------------------- 
     580       ! X) remove water from the snow 
     581       !------------------------------- 
     582       !------------------------------------------------------------------------ 
     583       ! total melt pond volume in category DOes not include snow volume 
     584       ! snow in melt ponds is not melted 
     585       !------------------------------------------------------------------------ 
     586  
     587       ! Calculate pond volume for lower categories 
     588       DO n=1,m_index-1 
     589          zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 
     590                   - (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n)) 
     591       END DO 
     592  
     593       ! Calculate pond volume for highest category = remaining pond volume 
     594 
     595       ! The following is completely unclear to Martin at least 
     596       ! Could we redefine properly and recode in a more readable way ? 
     597 
     598       ! m_index = last category with melt pond 
     599 
     600       IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 
     601 
     602       IF (m_index > 1) THEN 
     603         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 
     604           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) !  
     605         ELSE 
     606           zvolpn(m_index) = z0 
     607           zhpondn(m_index) = z0 
     608           zapondn(m_index) = z0 
     609           ! If remaining pond volume is negative reduce pond volume of  
     610           ! lower category 
     611           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &  
     612             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 
     613                                + zvolp 
     614         END IF 
     615       END IF 
     616  
     617       DO n=1,m_index 
     618          IF (zapondn(n) > epsi10) THEN 
     619              zhpondn(n) = zvolpn(n) / zapondn(n) 
     620          ELSE 
     621             dvolp = dvolp + zvolpn(n) 
     622             zhpondn(n) = z0 
     623             zvolpn(n) = z0 
     624             zapondn(n) = z0 
     625          end IF 
     626       END DO 
     627       DO n = m_index+1, jpl 
     628          zhpondn(n) = z0 
     629          zapondn(n) = z0 
     630          zvolpn (n) = z0 
     631       END DO 
     632  
     633    END SUBROUTINE lim_mp_area 
     634 
     635!OLI_CODE    
     636!OLI_CODE  
     637!OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, & 
     638!OLI_CODE                          zvolp, cum_max_vol,                & 
     639!OLI_CODE                          hpond, m_index) 
     640!OLI_CODE       !!------------------------------------------------------------------- 
     641!OLI_CODE       !!                ***  ROUTINE calc_hpond  *** 
     642!OLI_CODE       !! 
     643!OLI_CODE       !! ** Purpose :   Compute melt pond depth 
     644!OLI_CODE       !!------------------------------------------------------------------- 
     645!OLI_CODE      
     646!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     647!OLI_CODE          aicen, & 
     648!OLI_CODE          asnon, & 
     649!OLI_CODE          hsnon, & 
     650!OLI_CODE          rhos,  & 
     651!OLI_CODE          alfan, & 
     652!OLI_CODE          cum_max_vol 
     653!OLI_CODE      
     654!OLI_CODE       REAL (wp), INTENT(IN) :: & 
     655!OLI_CODE          zvolp 
     656!OLI_CODE      
     657!OLI_CODE       REAL (wp), INTENT(OUT) :: & 
     658!OLI_CODE          hpond 
     659!OLI_CODE      
     660!OLI_CODE       INTEGER, INTENT(OUT) :: & 
     661!OLI_CODE          m_index 
     662!OLI_CODE      
     663!OLI_CODE       INTEGER :: n, ns 
     664!OLI_CODE      
     665!OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: & 
     666!OLI_CODE          hitl, & 
     667!OLI_CODE          aicetl 
     668!OLI_CODE      
     669!OLI_CODE       REAL (wp) :: & 
     670!OLI_CODE          rem_vol, & 
     671!OLI_CODE          area, & 
     672!OLI_CODE          vol, & 
     673!OLI_CODE          tmp, & 
     674!OLI_CODE          z0   = 0.0_wp,    &    
     675!OLI_CODE          epsi10 = 1.0e-11_wp 
     676!OLI_CODE      
     677!OLI_CODE       !---------------------------------------------------------------- 
     678!OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained?  
     679!OLI_CODE       !---------------------------------------------------------------- 
     680!OLI_CODE      
     681!OLI_CODE       IF (zvolp < epsi10) THEN 
     682!OLI_CODE        hpond = z0 
     683!OLI_CODE        m_index = 0 
     684!OLI_CODE       ELSE 
     685!OLI_CODE         
     686!OLI_CODE        !---------------------------------------------------------------- 
     687!OLI_CODE        ! Calculate the category where water fills up to  
     688!OLI_CODE        !---------------------------------------------------------------- 
     689!OLI_CODE         
     690!OLI_CODE        !----------| 
     691!OLI_CODE        !          | 
     692!OLI_CODE        !          | 
     693!OLI_CODE        !          |----------|                                     -- -- 
     694!OLI_CODE        !__________|__________|_________________________________________ ^ 
     695!OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled 
     696!OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     697!OLI_CODE        !          |          |          |              | 
     698!OLI_CODE        !          |          |          |              |hpond 
     699!OLI_CODE        !          |          |          |----------|   |     |------- 
     700!OLI_CODE        !          |          |          |          |   |     | 
     701!OLI_CODE        !          |          |          |          |---v-----| 
     702!OLI_CODE        !          |          | m_index  |          |         | 
     703!OLI_CODE        !------------------------------------------------------------- 
     704!OLI_CODE         
     705!OLI_CODE        m_index = 0  ! 1:m_index categories have water in them 
     706!OLI_CODE        DO n = 1, jpl 
     707!OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN 
     708!OLI_CODE              m_index = n 
     709!OLI_CODE              IF (n == 1) THEN 
     710!OLI_CODE                 rem_vol = zvolp 
     711!OLI_CODE              ELSE 
     712!OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1) 
     713!OLI_CODE              END IF 
     714!OLI_CODE              exit ! to break out of the loop 
     715!OLI_CODE           END IF 
     716!OLI_CODE        END DO 
     717!OLI_CODE        m_index = min(jpl-1, m_index) 
     718!OLI_CODE         
     719!OLI_CODE        !---------------------------------------------------------------- 
     720!OLI_CODE        ! semi-filled layer may have m_index different snow in it 
     721!OLI_CODE        !---------------------------------------------------------------- 
     722!OLI_CODE         
     723!OLI_CODE        !-----------------------------------------------------------  ^ 
     724!OLI_CODE        !                                                             |  alfan(m_index+1) 
     725!OLI_CODE        !                                                             | 
     726!OLI_CODE        !hitl(3)-->                             |----------|          | 
     727!OLI_CODE        !hitl(2)-->                |------------| * * * * *|          | 
     728!OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     729!OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^ 
     730!OLI_CODE        !                various snow from lower categories          |  |alfa(m_index) 
     731!OLI_CODE         
     732!OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories 
     733!OLI_CODE        ! aicetl - area of each snow depth in this layer 
     734!OLI_CODE         
     735!OLI_CODE        hitl(:) = z0 
     736!OLI_CODE        aicetl(:) = z0 
     737!OLI_CODE        DO n = 1, m_index 
     738!OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     739!OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0) 
     740!OLI_CODE           aicetl(n) = asnon(n) 
     741!OLI_CODE            
     742!OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     743!OLI_CODE        END DO 
     744!OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     745!OLI_CODE        aicetl(m_index+1) = z0 
     746!OLI_CODE         
     747!OLI_CODE        !---------------------------------------------------------------- 
     748!OLI_CODE        ! reorder array according to hitl  
     749!OLI_CODE        ! snow heights not necessarily in height order 
     750!OLI_CODE        !---------------------------------------------------------------- 
     751!OLI_CODE         
     752!OLI_CODE        DO ns = 1, m_index+1 
     753!OLI_CODE           DO n = 0, m_index - ns + 1 
     754!OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     755!OLI_CODE                 tmp = hitl(n) 
     756!OLI_CODE                 hitl(n) = hitl(n+1) 
     757!OLI_CODE                 hitl(n+1) = tmp 
     758!OLI_CODE                 tmp = aicetl(n) 
     759!OLI_CODE                 aicetl(n) = aicetl(n+1) 
     760!OLI_CODE                 aicetl(n+1) = tmp 
     761!OLI_CODE              END IF 
     762!OLI_CODE           END DO 
     763!OLI_CODE        END DO 
     764!OLI_CODE         
     765!OLI_CODE        !---------------------------------------------------------------- 
     766!OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     767!OLI_CODE        !---------------------------------------------------------------- 
     768!OLI_CODE         
     769!OLI_CODE        !hitl(3)---------------------------------------------------------------- 
     770!OLI_CODE        !                                                       | * * * * * * * *   
     771!OLI_CODE        !                                                       |* * * * * * * * *  
     772!OLI_CODE        !hitl(2)---------------------------------------------------------------- 
     773!OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * *   
     774!OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * *  
     775!OLI_CODE        !hitl(1)---------------------------------------------------------------- 
     776!OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *   
     777!OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *  
     778!OLI_CODE        !hitl(0)---------------------------------------------------------------- 
     779!OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)             
     780!OLI_CODE         
     781!OLI_CODE        ! move up over layers incrementing volume 
     782!OLI_CODE        DO n = 1, m_index+1 
     783!OLI_CODE            
     784!OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     785!OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     786!OLI_CODE            
     787!OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     788!OLI_CODE            
     789!OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     790!OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - & 
     791!OLI_CODE                           alfan(1) 
     792!OLI_CODE              exit 
     793!OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth 
     794!OLI_CODE              rem_vol = rem_vol - vol 
     795!OLI_CODE           END IF 
     796!OLI_CODE            
     797!OLI_CODE        END DO 
     798!OLI_CODE         
     799!OLI_CODE       END IF 
     800!OLI_CODE      
     801!OLI_CODE    END SUBROUTINE calc_hpond 
     802!OLI_CODE    
     803!OLI_CODE  
     804    SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm) 
     805       !!------------------------------------------------------------------- 
     806       !!                ***  ROUTINE lim_mp_perm *** 
     807       !! 
     808       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     809       !!                and its permeability 
     810       !!------------------------------------------------------------------- 
     811       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     812          ticen,  & ! energy of melting for each ice layer (J/m2) 
     813          salin 
     814  
     815       REAL (wp), INTENT(IN) :: & 
     816          vicen     ! ice volume 
     817      
     818       REAL (wp), INTENT(OUT) :: & 
     819          perm      ! permeability 
     820  
     821       REAL (wp) ::   & 
     822          Sbr        ! brine salinity 
     823  
     824       REAL (wp), DIMENSION(nlay_i) ::   & 
     825          Tin, &    ! ice temperature 
     826          phi       ! liquid fraction 
     827  
     828       INTEGER :: k 
     829      
     830       REAL (wp) :: & 
     831          c2    = 2.0_wp 
     832   
     833       !----------------------------------------------------------------- 
     834       ! Compute ice temperatures from enthalpies using quadratic formula 
     835       !----------------------------------------------------------------- 
     836  
     837       DO k = 1,nlay_i 
     838          Tin(k) = ticen(k)  
     839       END DO 
     840  
     841       !----------------------------------------------------------------- 
     842       ! brine salinity and liquid fraction 
     843       !----------------------------------------------------------------- 
     844  
     845       IF (maxval(Tin-rtt) <= -c2) THEN 
     846  
     847          DO k = 1,nlay_i 
     848             Sbr = - 1.2_wp                 & 
     849                   -21.8_wp     * (Tin(k)-rtt)    & 
     850                   - 0.919_wp   * (Tin(k)-rtt)**2 & 
     851                   - 0.01878_wp * (Tin(k)-rtt)**3 
     852             phi(k) = salin(k)/Sbr ! liquid fraction 
     853          END DO ! k 
     854         
     855       ELSE 
     856  
     857          DO k = 1,nlay_i 
     858             Sbr = -17.6_wp    * (Tin(k)-rtt)    & 
     859                   - 0.389_wp  * (Tin(k)-rtt)**2 & 
     860                   - 0.00362_wp* (Tin(k)-rtt)**3 
     861             phi(k) = salin(k)/Sbr ! liquid fraction 
     862          END DO 
     863  
     864       END IF 
     865      
     866       !----------------------------------------------------------------- 
     867       ! permeability 
     868       !----------------------------------------------------------------- 
     869  
     870       perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking 
     871                                            ! bastard of Golden) 
     872      
     873    END SUBROUTINE lim_mp_perm 
     874!OLI_CODE    
     875!OLI_CODE #else 
     876!OLI_CODE    !!---------------------------------------------------------------------- 
     877!OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model 
     878!OLI_CODE    !!---------------------------------------------------------------------- 
     879!OLI_CODE CONTAINS 
     880!OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine 
     881!OLI_CODE    END SUBROUTINE lim_mp_init    
     882!OLI_CODE    SUBROUTINE lim_mp               ! Empty routine 
     883!OLI_CODE    END SUBROUTINE lim_mp 
     884!OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine 
     885!OLI_CODE    END SUBROUTINE compute_mp_topo         
     886!OLI_CODE    SUBROUTINE pond_area            ! Empty routine 
     887!OLI_CODE    END SUBROUTINE pond_area    
     888!OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine 
     889!OLI_CODE    END SUBROUTINE calc_hpond    
     890!OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine 
     891!OLI_CODE    END SUBROUTINE permeability_phy    
     892!OLI_CODE #endif 
     893!OLI_CODE    !!====================================================================== 
     894!OLI_CODE END MODULE limmp_topo 
     895 
    147896 
    148897!OLI_CODE MODULE limmp_topo 
     
    3651114!OLI_CODE          z25      = 0.25_wp,  & 
    3661115!OLI_CODE          z5       = 0.5_wp,   & 
    367 !OLI_CODE          zuny     = 1.0e-11_wp 
     1116!OLI_CODE          epsi10     = 1.0e-11_wp 
    3681117!OLI_CODE       !--------------------------------------------------------------- 
    3691118!OLI_CODE       ! initialize 
     
    4091158!OLI_CODE       DO i = 1, jpi 
    4101159!OLI_CODE          zhi = z0 
    411 !OLI_CODE          IF (aice(i,j) > zuny) zhi = vice(i,j)/aice(i,j) 
     1160!OLI_CODE          IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j) 
    4121161!OLI_CODE          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 
    4131162!OLI_CODE             zvolp(i,j) > zmin_volp*aice(i,j)) THEN 
     
    4481197!OLI_CODE          DO n = 1, jpl-1 
    4491198!OLI_CODE                         
    450 !OLI_CODE             IF (zvuin(i,j,n) > zuny) THEN 
     1199!OLI_CODE             IF (zvuin(i,j,n) > epsi10) THEN 
    4511200!OLI_CODE  
    4521201!OLI_CODE          !---------------------------------------------------------------- 
     
    4571206!OLI_CODE  
    4581207!OLI_CODE                   zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n)) 
    459 !OLI_CODE                   IF (zdvice > zuny) THEN 
     1208!OLI_CODE                   IF (zdvice > epsi10) THEN 
    4601209!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) - zdvice 
    4611210!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice 
     
    4631212!OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   + zdvice 
    4641213!OLI_CODE                         
    465 !OLI_CODE                      IF (zvuin(i,j,n) < zuny .and. zvolpn(i,j,n) > puny) THEN 
     1214!OLI_CODE                      IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN 
    4661215!OLI_CODE                         ! ice lid melted and category is pond covered 
    4671216!OLI_CODE                         zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n) 
     
    4751224!OLI_CODE          ! freezing: existing upper ice layer grows 
    4761225!OLI_CODE          !---------------------------------------------------------------- 
    477 !OLI_CODE                ELSE IF (zvolpn(i,j,n) > zuny) THEN ! zTavg <= zTp 
     1226!OLI_CODE                ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp 
    4781227!OLI_CODE  
    4791228!OLI_CODE                 ! dIFferential growth of base of surface floating ice layer 
     
    4841233!OLI_CODE  
    4851234!OLI_CODE                   zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))    
    486 !OLI_CODE                   IF (zdvice > zuny) THEN 
     1235!OLI_CODE                   IF (zdvice > epsi10) THEN 
    4871236!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) + zdvice 
    4881237!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 
     
    4981247!OLI_CODE          ! note: albedo does not change 
    4991248!OLI_CODE          !---------------------------------------------------------------- 
    500 !OLI_CODE             ELSE ! zvuin < zuny 
     1249!OLI_CODE             ELSE ! zvuin < epsi10 
    5011250!OLI_CODE                      
    5021251!OLI_CODE                ! thickness of newly formed ice 
     
    5061255!OLI_CODE                zdHui = max(-fsurf(i,j)*rdt_ice*zr1_rlfus, z0) 
    5071256!OLI_CODE                zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))   
    508 !OLI_CODE                IF (zdvice > zuny) THEN 
     1257!OLI_CODE                IF (zdvice > epsi10) THEN 
    5091258!OLI_CODE                   zvuin (i,j,n) = zdvice 
    5101259!OLI_CODE                   zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 
     
    5281277!OLI_CODE       DO i = 1, jpi 
    5291278!OLI_CODE          DO n = 1, jpl 
    530 !OLI_CODE             IF (aicen(i,j,n) > zuny .and. zvolpn(i,j,n) < puny & 
    531 !OLI_CODE                                     .and. zvuin (i,j,n) > zuny) THEN 
     1279!OLI_CODE             IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny & 
     1280!OLI_CODE                                     .and. zvuin (i,j,n) > epsi10) THEN 
    5321281!OLI_CODE                kcells(n) = kcells(n) + 1 
    5331282!OLI_CODE                indxij    = kcells(n) 
     
    5531302!OLI_CODE          DO j = 1, jpj 
    5541303!OLI_CODE             DO i = 1, jpi 
    555 !OLI_CODE                IF (zapondn(i,j,n) > zuny) THEN 
     1304!OLI_CODE                IF (zapondn(i,j,n) > epsi10) THEN 
    5561305!OLI_CODE                   h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n) 
    5571306!OLI_CODE                ELSE 
     
    5591308!OLI_CODE                   h_il(i,j,n) = z0 
    5601309!OLI_CODE                END IF 
    561 !OLI_CODE                IF (aicen(i,j,n) > zuny) THEN 
     1310!OLI_CODE                IF (aicen(i,j,n) > epsi10) THEN 
    5621311!OLI_CODE                   a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * & 
    5631312!OLI_CODE                         (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n)))) 
     
    6371386!OLI_CODE          p4        = 0.4_wp    , & 
    6381387!OLI_CODE          p6        = 0.6_wp    , & 
    639 !OLI_CODE          zuny      = 1.0e-11_wp 
     1388!OLI_CODE          epsi10      = 1.0e-11_wp 
    6401389!OLI_CODE           
    6411390!OLI_CODE      !-----------| 
     
    6711420!OLI_CODE          zhpondn(n) = z0 
    6721421!OLI_CODE  
    673 !OLI_CODE          IF (aicen(n) < zuny)  THEN 
     1422!OLI_CODE          IF (aicen(n) < epsi10)  THEN 
    6741423!OLI_CODE             hicen(n) =  z0  
    6751424!OLI_CODE             hsnon(n) = z0 
     
    7421491!OLI_CODE       !---------------------------------------------------------------- 
    7431492!OLI_CODE       IF (zvolp >= cum_max_vol(jpl)) THEN 
    744 !OLI_CODE          drain = zvolp - cum_max_vol(jpl) + zuny 
     1493!OLI_CODE          drain = zvolp - cum_max_vol(jpl) + epsi10 
    7451494!OLI_CODE          zvolp = zvolp - drain 
    7461495!OLI_CODE          dvolp = drain 
    747 !OLI_CODE          IF (zvolp < zuny) THEN 
     1496!OLI_CODE          IF (zvolp < epsi10) THEN 
    7481497!OLI_CODE             dvolp = dvolp + zvolp 
    7491498!OLI_CODE             zvolp = z0 
     
    7891538!OLI_CODE             dvolp = dvolp + min(drain, zvolp) 
    7901539!OLI_CODE             zvolp = max(zvolp - drain, z0) 
    791 !OLI_CODE             IF (zvolp < zuny) THEN 
     1540!OLI_CODE             IF (zvolp < epsi10) THEN 
    7921541!OLI_CODE                dvolp = dvolp + zvolp 
    7931542!OLI_CODE                zvolp = z0 
     
    8311580!OLI_CODE           ! If remaining pond volume is negative reduce pond volume of  
    8321581!OLI_CODE           ! lower category 
    833 !OLI_CODE           IF (zvolp+zuny < sum(zvolpn(1:m_index-1))) &  
     1582!OLI_CODE           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &  
    8341583!OLI_CODE             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 
    8351584!OLI_CODE                                + zvolp 
     
    8381587!OLI_CODE  
    8391588!OLI_CODE       DO n=1,m_index 
    840 !OLI_CODE          IF (zapondn(n) > zuny) THEN 
     1589!OLI_CODE          IF (zapondn(n) > epsi10) THEN 
    8411590!OLI_CODE              zhpondn(n) = zvolpn(n) / zapondn(n) 
    8421591!OLI_CODE          ELSE 
     
    8941643!OLI_CODE          tmp, & 
    8951644!OLI_CODE          z0   = 0.0_wp,    &    
    896 !OLI_CODE          zuny = 1.0e-11_wp 
     1645!OLI_CODE          epsi10 = 1.0e-11_wp 
    8971646!OLI_CODE      
    8981647!OLI_CODE       !---------------------------------------------------------------- 
     
    9001649!OLI_CODE       !---------------------------------------------------------------- 
    9011650!OLI_CODE      
    902 !OLI_CODE       IF (zvolp < zuny) THEN 
     1651!OLI_CODE       IF (zvolp < epsi10) THEN 
    9031652!OLI_CODE        hpond = z0 
    9041653!OLI_CODE        m_index = 0 
     
    11141863!OLI_CODE END MODULE limmp_topo 
    11151864!OLI_CODE  
     1865#else 
     1866   !!---------------------------------------------------------------------- 
     1867   !!   Default option          Empty module           NO LIM sea-ice model 
     1868   !!---------------------------------------------------------------------- 
     1869CONTAINS 
     1870   SUBROUTINE lim_mp_init     ! Empty routine 
     1871   END SUBROUTINE lim_mp_init 
     1872   SUBROUTINE lim_mp          ! Empty routine 
     1873   END SUBROUTINE lim_mp      
     1874   SUBROUTINE lim_mp_topo     ! Empty routine 
     1875   END SUBROUTINE lim_mp_topo 
     1876   SUBROUTINE lim_mp_area     ! Empty routine 
     1877   END SUBROUTINE lim_mp_area 
     1878   SUBROUTINE lim_mp_perm     ! Empty routine 
     1879   END SUBROUTINE lim_mp_perm 
     1880#endif  
     1881 
     1882   !!====================================================================== 
     1883END MODULE limmp  
  • branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r7319 r7325  
    261261         ! --- zap this if no ice thermo --- ! 
    262262         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics       
     263 
     264         ! MV MP 2016 
     265         IF ( ln_limMP)               CALL lim_mp( kt )         ! -- Melt ponds 
     266         ! END MV MP 2016 
     267 
    263268         IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections 
    264269         ! --- 
Note: See TracChangeset for help on using the changeset viewer.