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.
limmp.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90 @ 8233

Last change on this file since 8233 was 8233, checked in by clem, 7 years ago

merge with dev_r6859_LIM3_meltponds@r8179

File size: 93.5 KB
Line 
1MODULE limmp 
2   !!======================================================================
3   !!                     ***  MODULE  limmp   ***
4   !!   Melt ponds
5   !!======================================================================
6   !! history :       ! Original code by Daniela Flocco and Adrian Turner
7   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
8   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
9   !!                 ! NB: Only lim_mp_cstt and lim_mp_cesm work in this
10   !!                       version
11   !!----------------------------------------------------------------------
12#if defined key_lim3
13   !!----------------------------------------------------------------------
14   !!   'key_lim3' :                                 LIM3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   lim_mp_init      : some initialization and namelist read
17   !!   lim_mp           : main calling routine
18   !!   lim_mp_topo      : main melt pond routine for the "topographic" formulation (FloccoFeltham)
19   !!   lim_mp_area      : ??? compute melt pond fraction per category
20   !!   lim_mp_perm      : computes permeability (should be a FUNCTION!)
21   !!   calc_hpond       : computes melt pond depth
22   !!   permeability_phy : computes permeability
23
24   !!----------------------------------------------------------------------
25   USE phycst           ! physical constants
26   USE dom_oce          ! ocean space and time domain
27!  USE sbc_ice          ! Surface boundary condition: ice   fields
28   USE ice              ! LIM-3 variables
29   USE lbclnk           ! lateral boundary conditions - MPP exchanges
30   USE lib_mpp          ! MPP library
31   USE wrk_nemo         ! work arrays
32   USE in_out_manager   ! I/O manager
33   USE lib_fortran      ! glob_sum
34   USE timing           ! Timing
35!  USE limcons          ! conservation tests
36!  USE limctl           ! control prints
37!  USE limvar
38
39!OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice
40!OLI_CODE    USE phycst
41!OLI_CODE    USE dom_ice
42!OLI_CODE    USE dom_oce
43!OLI_CODE    USE sbc_oce
44!OLI_CODE    USE sbc_ice
45!OLI_CODE    USE par_ice
46!OLI_CODE    USE par_oce
47!OLI_CODE    USE ice
48!OLI_CODE    USE thd_ice
49!OLI_CODE    USE in_out_manager
50!OLI_CODE    USE lbclnk
51!OLI_CODE    USE lib_mpp
52!OLI_CODE
53!OLI_CODE    IMPLICIT NONE
54!OLI_CODE    PRIVATE
55!OLI_CODE
56!OLI_CODE    PUBLIC   lim_mp_init
57!OLI_CODE    PUBLIC   lim_mp
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   lim_mp_init    ! routine called by sbcice_lim.F90
63   PUBLIC   lim_mp         ! routine called by sbcice_lim.F90
64
65   !! * Substitutions
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
69   !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE lim_mp_init 
75      !!-------------------------------------------------------------------
76      !!                  ***  ROUTINE lim_mp_init   ***
77      !!
78      !! ** Purpose : Physical constants and parameters linked to melt ponds
79      !!      over sea ice
80      !!
81      !! ** Method  :  Read the namicemp  namelist and check the melt pond 
82      !!       parameter values called at the first timestep (nit000)
83      !!
84      !! ** input   :   Namelist namicemp 
85      !!-------------------------------------------------------------------
86      INTEGER  ::   ios                 ! Local integer output status for namelist read
87      NAMELIST/namicemp/  ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd
88      !!-------------------------------------------------------------------
89
90      REWIND( numnam_ice_ref )              ! Namelist namicemp  in reference namelist : Melt Ponds 
91      READ  ( numnam_ice_ref, namicemp, IOSTAT = ios, ERR = 901)
92901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp  in reference namelist', lwp )
93
94      REWIND( numnam_ice_cfg )              ! Namelist namicemp  in configuration namelist : Melt Ponds
95      READ  ( numnam_ice_cfg, namicemp, IOSTAT = ios, ERR = 902 )
96902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp in configuration namelist', lwp )
97      IF(lwm) WRITE ( numoni, namicemp )
98     
99      IF(lwp) THEN                        ! control print
100         WRITE(numout,*)
101         WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds'
102         WRITE(numout,*) '~~~~~~~~~~~~'
103         WRITE(numout,*) '   Active melt ponds                                           ln_pnd        = ', ln_pnd
104         WRITE(numout,*) '   Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad
105         WRITE(numout,*) '   Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
106         WRITE(numout,*) '   Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme
107         WRITE(numout,*) '   Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
108         WRITE(numout,*) '   Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd
109      ENDIF
110
111      IF ( .NOT. ln_pnd ) THEN
112         WRITE(numout,*)
113         WRITE(numout,*) ' Melt ponds are not activated '
114         WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '
115         WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '
116         ln_pnd_rad    = .FALSE.
117         ln_pnd_fw     = .FALSE.
118         nn_pnd_scheme = 0
119         rn_apnd       = 0._wp
120         rn_hpnd       = 0._wp
121
122         IF(lwp) THEN                     ! control print
123            WRITE(numout,*) '    Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad
124            WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
125            WRITE(numout,*) '    Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme
126            WRITE(numout,*) '    Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
127            WRITE(numout,*) '    Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd
128         ENDIF
129      ENDIF
130
131      IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN
132         WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false ' 
133         ln_pnd_fw     = .FALSE.
134         IF(lwp) THEN                     ! control print
135            WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
136         ENDIF
137      ENDIF
138
139      IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN
140         WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
141         WRITE(numout,*) ' Run aborted '
142         CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' )
143      ENDIF
144
145      !
146   END SUBROUTINE lim_mp_init
147   
148
149   SUBROUTINE lim_mp( kt )
150      !!-------------------------------------------------------------------
151      !!               ***  ROUTINE lim_mp   ***
152      !!               
153      !! ** Purpose :   change melt pond fraction
154      !!               
155      !! ** Method  :   brutal force
156      !!
157      !! ** Action  : -
158      !!              -
159      !!------------------------------------------------------------------------------------
160
161      INTEGER, INTENT(in) :: kt    ! number of iteration
162      INTEGER  ::   ji, jj, jl     ! dummy loop indices
163
164      !!-------------------------------------------------------------------
165
166      IF( nn_timing == 1 )  CALL timing_start('lim_mp')
167
168      SELECT CASE ( nn_pnd_scheme )
169
170         CASE (0)
171
172            CALL lim_mp_cstt       ! staircase melt ponds
173
174         CASE (1)
175
176            CALL lim_mp_cesm       ! empirical melt ponds
177
178         CASE (2)
179
180            CALL lim_mp_topo   &   ! topographic melt ponds
181                      &          (at_i, a_i,                                       &
182                      &          vt_i, v_i, v_s,            t_i, s_i, a_ip_frac,   &
183                      &          h_ip,     t_su)
184
185      END SELECT
186
187      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
188
189   END SUBROUTINE lim_mp 
190
191   SUBROUTINE lim_mp_cstt 
192       !!-------------------------------------------------------------------
193       !!                ***  ROUTINE lim_mp_cstt  ***
194       !!
195       !! ** Purpose    : Compute melt pond evolution
196       !!
197       !! ** Method     : Melt pond fraction and thickness are prescribed
198       !!                 to non-zero values when t_su = 0C
199       !!
200       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
201       !!               
202       !! ** Note       : Coupling with such melt ponds is only radiative
203       !!                 Advection, ridging, rafting... are bypassed
204       !!
205       !! ** References : Bush, G.W., and Trump, D.J. (2017)
206       !!   
207       !!-------------------------------------------------------------------
208       INTEGER                             :: ji, jj, jl
209       REAL(wp)                            :: z1_jpl            ! 1/jpl
210       !!-------------------------------------------------------------------
211
212       z1_jpl     = 1. / FLOAT(jpl)
213
214       WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
215          a_ip_frac = rn_apnd
216          h_ip      = rn_hpnd   
217          v_ip      = a_ip_frac * a_i * h_ip 
218          a_ip      = a_ip_frac * a_i
219       ELSE WHERE
220          a_ip      = 0._wp
221          h_ip      = 0._wp
222          v_ip      = 0._wp
223          a_ip_frac = 0._wp
224       END WHERE
225
226       wfx_pnd(:,:) = 0._wp
227
228   END SUBROUTINE lim_mp_cstt
229
230   SUBROUTINE lim_mp_cesm
231       !!-------------------------------------------------------------------
232       !!                ***  ROUTINE lim_mp_cesm  ***
233       !!
234       !! ** Purpose    : Compute melt pond evolution
235       !!
236       !! ** Method     : Empirical method. A fraction of meltwater is accumulated
237       !!                 in pond volume. It is then released exponentially when
238       !!                 surface is freezing.
239       !!
240       !! ** Tunable parameters : (no real expertise yet, ideas?)
241       !!
242       !! ** Note       : Stolen from CICE for quick test of the melt pond
243       !!                 radiation and freshwater interfaces
244       !!                 Coupling can be radiative AND freshwater
245       !!                 Advection, ridging, rafting are called
246       !!
247       !! ** References : Holland, M. M. et al (J Clim 2012)
248       !!   
249       !!-------------------------------------------------------------------
250
251       INTEGER, POINTER, DIMENSION(:)      :: indxi             ! compressed indices for cells with ice melting
252       INTEGER, POINTER, DIMENSION(:)      :: indxj             !
253
254       REAL(wp), POINTER, DIMENSION(:,:)   :: zwfx_mlw          ! available meltwater for melt ponding
255       REAL(wp), POINTER, DIMENSION(:,:,:) :: zrfrac            ! fraction of available meltwater retained for melt ponding
256
257       REAL(wp), PARAMETER                 :: zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
258       REAL(wp), PARAMETER                 :: zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
259       REAL(wp), PARAMETER                 :: zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
260       REAL(wp), PARAMETER                 :: zpnd_aspect = 0.8_wp ! pond aspect ratio
261
262       REAL(wp)                            :: zhi               ! dummy ice thickness
263       REAL(wp)                            :: zhs               ! dummy snow depth
264       REAL(wp)                            :: zTp               ! reference temperature
265       REAL(wp)                            :: zdTs              ! dummy temperature difference
266       REAL(wp)                            :: z1_rhofw          ! inverse freshwater density
267       REAL(wp)                            :: z1_zpnd_aspect    ! inverse pond aspect ratio
268       REAL(wp)                            :: zvpold            ! dummy pond volume
269
270       INTEGER                             :: ji, jj, jl, ij    ! loop indices
271       INTEGER                             :: icells            ! size of dummy array
272
273       !!-------------------------------------------------------------------
274
275        CALL wrk_alloc( jpi*jpj, indxi, indxj)
276        CALL wrk_alloc( jpi,jpj,     zwfx_mlw )
277        CALL wrk_alloc( jpi,jpj,jpl, zrfrac   )
278
279        z1_rhofw       = 1. / rhofw 
280        z1_zpnd_aspect = 1. / zpnd_aspect
281        zTp            = -2. 
282
283        a_ip_frac(:,:,:) = 0._wp
284        h_ip     (:,:,:) = 0._wp
285 
286        !------------------------------------------------------------------
287        ! Available melt water for melt ponding and corresponding fraction
288        !------------------------------------------------------------------
289 
290        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
291
292                ! NB: zwfx_mlw can be slightly negative for very small values (why?)
293                ! This can in some occasions give negative
294                ! v_ip in the first category, which then gives crazy pond
295                ! fractions and crashes the code as soon as the melt-pond
296                ! radiative coupling is activated
297                ! if we understand and remove why wfx_sum or wfx_snw could be
298                ! negative, then, we can remove the MAX
299                ! NB: I now changed to wfx_snw_sum, this may fix the problem.
300                ! We should check
301 
302        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:) 
303 
304        DO jl = 1, jpl   
305
306           !------------------------------------------------------------------------------
307           ! Identify grid cells where ponds should be updated (can probably be improved)
308           !------------------------------------------------------------------------------
309 
310           indxi(:) = 0
311           indxj(:) = 0
312           icells   = 0
313 
314           DO jj = 1, jpj
315             DO ji = 1, jpi
316                IF ( a_i(ji,jj,jl) > epsi10 ) THEN
317                   icells = icells + 1
318                   indxi(icells) = ji
319                   indxj(icells) = jj
320                ENDIF
321             END DO                 ! ji
322          END DO                    ! jj
323 
324          DO ij = 1, icells
325 
326             ji = indxi(ij)
327             jj = indxj(ij)
328 
329             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
330             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
331 
332             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
333 
334                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
335                v_ip(ji,jj,jl)      = 0._wp
336                a_ip_frac(ji,jj,jl) = 0._wp
337                h_ip(ji,jj,jl)      = 0._wp
338 
339                IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
340                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
341 
342 
343             ELSE                        !--- Update pond characteristics
344 
345                !--- Add retained melt water to melt ponds
346                ! v_ip should never be positive, otherwise code crashes
347                ! MV: as far as I saw, UM5 can create very small negative v_ip values
348                ! hence I added the max, which was not required with Prather (1 yr run)
349                v_ip(ji,jj,jl)      = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice
350 
351                !--- Shrink pond due to refreezing
352                zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
353               
354                zvpold              = v_ip(ji,jj,jl)
355 
356                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp )
357 
358                !--- Dump meltwater due to refreezing ( of course this is wrong
359                !--- but this parameterization is too simple )
360                IF ( ln_pnd_fw ) &
361                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
362 
363                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
364                   !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
365 
366                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
367 
368                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
369 
370             !-----------------------------------------------------------
371             ! Limit pond depth
372             !-----------------------------------------------------------
373             ! The original version has pond depth limitation, which I did not
374             ! keep here. Maybe we need it later on
375             !
376             
377             ENDIF
378 
379           END DO
380 
381        END DO ! jpl
382 
383        !--- Remove retained meltwater from surface fluxes
384 
385        IF ( ln_pnd_fw ) THEN
386 
387            wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
388
389            wfx_sum(:,:)     = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
390
391        ENDIF
392
393       CALL wrk_dealloc( jpi*jpj, indxi, indxj)
394       CALL wrk_dealloc( jpi,jpj,     zwfx_mlw )
395       CALL wrk_dealloc( jpi,jpj,jpl, zrfrac   )
396
397   END SUBROUTINE lim_mp_cesm
398
399   SUBROUTINE lim_mp_topo    (aice,      aicen,     &
400                              vice,      vicen,     &
401                              vsnon,                &
402                              ticen,     salin,     &
403                              a_ip_frac, h_ip,      &
404                                         Tsfc )
405       !!-------------------------------------------------------------------
406       !!                ***  ROUTINE lim_mp_topo  ***
407       !!
408       !! ** Purpose :   Compute melt pond evolution based on the ice
409       !!                topography as inferred from the ice thickness
410       !!                distribution. 
411       !!
412       !! ** Method  :   This code is initially based on Flocco and Feltham
413       !!                (2007) and Flocco et al. (2010). More to come...
414       !!
415       !! ** Tunable parameters :
416       !!
417       !! ** Note :
418       !!
419       !! ** References
420       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond
421       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:
422       !!    10.1029/2006JC003836.
423       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of
424       !!    a physically based melt pond scheme into the sea ice component of a
425       !!    climate model.  J. Geophys. Res. 115, C08012,
426       !!    doi: 10.1029/2009JC005568.
427       !!   
428       !!-------------------------------------------------------------------
429 
430       REAL (wp), DIMENSION (jpi,jpj), &
431          INTENT(IN) :: &
432          aice, &    ! total ice area fraction
433          vice       ! total ice volume (m)
434 
435       REAL (wp), DIMENSION (jpi,jpj,jpl), &
436          INTENT(IN) :: &
437          aicen, &   ! ice area fraction, per category
438          vsnon, &   ! snow volume, per category (m)
439          vicen      ! ice volume, per category (m)
440 
441       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
442          INTENT(IN) :: &
443          ticen, &   ! ice enthalpy, per category
444          salin
445 
446       REAL (wp), DIMENSION (jpi,jpj,jpl), &
447          INTENT(INOUT) :: &
448          a_ip_frac , &   ! pond area fraction of ice, per ice category
449          h_ip       ! pond depth, per ice category (m)
450 
451       REAL (wp), DIMENSION (jpi,jpj,jpl), &
452          INTENT(IN) :: &
453          Tsfc       ! snow/sea ice surface temperature
454 
455       ! local variables
456       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
457          zTsfcn, & ! ice/snow surface temperature (C)
458          zvolpn, & ! pond volume per unit area, per category (m)
459          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
460 
461       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
462          zapondn,& ! pond area fraction, per category
463          zhpondn   ! pond depth, per category (m)
464 
465       REAL (wp), DIMENSION (jpi,jpj) :: &
466          zvolp       ! total volume of pond, per unit area of pond (m)
467 
468       REAL (wp) :: &
469          zhi,    & ! ice thickness (m)
470          zdHui,  & ! change in thickness of ice lid (m)
471          zomega, & ! conduction
472          zdTice, & ! temperature difference across ice lid (C)
473          zdvice, & ! change in ice volume (m)
474          zTavg,  & ! mean surface temperature across categories (C)
475          zTp,    & ! pond freezing temperature (C)
476          zdvn      ! change in melt pond volume for fresh water budget
477       INTEGER, DIMENSION (jpi*jpj) :: &
478          indxi, indxj    ! compressed indices for cells with ice melting
479 
480       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
481 
482       INTEGER, DIMENSION (jpl) :: &
483          kcells          ! cells where ice lid combines with vice
484 
485       INTEGER, DIMENSION (jpi*jpj,jpl) :: &
486          indxii, indxjj  ! i,j indices for kcells loop
487 
488       REAL (wp), parameter :: &
489          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m)
490          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C)
491          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
492          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
493          z0       = 0._wp,    & 
494          zTimelt   = 0._wp,    &
495          z01      = 0.01_wp,  &
496          z25      = 0.25_wp,  &
497          z5       = 0.5_wp
498
499       !---------------------------------------------------------------
500       ! Initialization
501       !---------------------------------------------------------------
502 
503       zhpondn(:,:,:) = 0._wp
504       zapondn(:,:,:) = 0._wp
505       indxii(:,:) = 0
506       indxjj(:,:) = 0
507       kcells(:)   = 0
508
509       zvolp(:,:) = wfx_sum(:,:) + wfx_snw_sum(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds
510       zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0                   ! Convert in Celsius
511 
512       ! The freezing temperature for meltponds is assumed slightly below 0C,
513       ! as if meltponds had a little salt in them.  The salt budget is not
514       ! altered for meltponds, but if it were then an actual pond freezing
515       ! temperature could be computed.
516 
517       ! zTp = zTimelt - zTd  ---> for lids
518 
519       !-----------------------------------------------------------------
520       ! Identify grid cells with ponds
521       !-----------------------------------------------------------------
522 
523       icells = 0
524       DO j = 1, jpj
525       DO i = 1, jpi
526          zhi = z0
527          IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j)
528          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
529             zvolp(i,j) > zmin_volp*aice(i,j)) THEN
530             icells = icells + 1
531             indxi(icells) = i
532             indxj(icells) = j
533          ELSE  ! remove ponds on thin ice
534             !fpond(i,j) = fpond(i,j) - zvolp(i,j)
535             zvolpn(i,j,:) = z0
536             zvuin (i,j,:) = z0
537             zvolp (i,j) = z0
538          END IF
539       END DO                     ! i
540       END DO                     ! j
541 
542       DO ij = 1, icells
543          i = indxi(ij)
544          j = indxj(ij)
545 
546          !--------------------------------------------------------------
547          ! calculate pond area and depth
548          !--------------------------------------------------------------
549          CALL lim_mp_area(aice(i,j),vice(i,j), &
550                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
551                    ticen(i,j,:,:), salin(i,j,:,:), &
552                    zvolpn(i,j,:), zvolp(i,j), &
553                    zapondn(i,j,:),zhpondn(i,j,:), zdvn)
554         ! outputs are
555         ! - zdvn
556         ! - zvolpn
557         ! - zvolp
558         ! - zapondn
559         ! - zhpondn
560
561          wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean
562 
563          ! mean surface temperature MV - why do we need that ? --> for the lid
564
565          ! zTavg = z0
566          ! DO n = 1, jpl
567          !   zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
568          ! END DO
569          ! zTavg = zTavg / aice(i,j)
570 
571       END DO ! ij
572 
573       !---------------------------------------------------------------
574       ! Update pond volume and fraction
575       !---------------------------------------------------------------
576 
577       a_ip(:,:,:) = zapondn(:,:,:)
578       v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:)
579       a_ip_frac(:,:,:) = 0._wp
580       h_ip     (:,:,:) = 0._wp
581
582    END SUBROUTINE lim_mp_topo
583
584    SUBROUTINE lim_mp_area(aice,vice, &
585                         aicen, vicen, vsnon, ticen, &
586                         salin, zvolpn, zvolp,         &
587                         zapondn,zhpondn,dvolp)
588
589       !!-------------------------------------------------------------------
590       !!                ***  ROUTINE lim_mp_area ***
591       !!
592       !! ** Purpose : Given the total volume of meltwater, update
593       !!              pond fraction (a_ip) and depth (should be volume)
594       !!
595       !! **
596       !!             
597       !!------------------------------------------------------------------
598
599       REAL (wp), INTENT(IN) :: &
600          aice,vice
601 
602       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
603          aicen, vicen, vsnon
604 
605       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
606          ticen, salin
607 
608       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
609          zvolpn
610 
611       REAL (wp), INTENT(INOUT) :: &
612          zvolp, dvolp
613 
614       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
615          zapondn, zhpondn
616 
617       INTEGER :: &
618          n, ns,   &
619          m_index, &
620          permflag
621 
622       REAL (wp), DIMENSION(jpl) :: &
623          hicen, &
624          hsnon, &
625          asnon, &
626          alfan, &
627          betan, &
628          cum_max_vol, &
629          reduced_aicen       
630 
631       REAL (wp), DIMENSION(0:jpl) :: &
632          cum_max_vol_tmp
633 
634       REAL (wp) :: &
635          hpond, &
636          drain, &
637          floe_weight, &
638          pressure_head, &
639          hsl_rel, &
640          deltah, &
641          perm, &
642          msno
643 
644       REAL (wp), parameter :: & 
645          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s
646          z0        = 0.0_wp    , &
647          c1        = 1.0_wp    , &
648          p4        = 0.4_wp    , &
649          p6        = 0.6_wp    , &
650          epsi10      = 1.0e-11_wp
651         
652      !-----------|
653      !           |
654      !           |-----------|
655      !___________|___________|______________________________________sea-level
656      !           |           |
657      !           |           |---^--------|
658      !           |           |   |        |
659      !           |           |   |        |-----------|              |-------
660      !           |           |   |alfan(n)|           |              |
661      !           |           |   |        |           |--------------|
662      !           |           |   |        |           |              |
663      !---------------------------v-------------------------------------------
664      !           |           |   ^        |           |              |
665      !           |           |   |        |           |--------------|
666      !           |           |   |betan(n)|           |              |
667      !           |           |   |        |-----------|              |-------
668      !           |           |   |        |
669      !           |           |---v------- |
670      !           |           |
671      !           |-----------|
672      !           |
673      !-----------|
674     
675       !-------------------------------------------------------------------
676       ! initialize
677       !-------------------------------------------------------------------
678 
679       DO n = 1, jpl
680 
681          zapondn(n) = z0
682          zhpondn(n) = z0
683 
684          !----------------------------------------
685          ! X) compute the effective snow fraction
686          !----------------------------------------
687          IF (aicen(n) < epsi10)  THEN
688             hicen(n) =  z0 
689             hsnon(n) = z0
690             reduced_aicen(n) = z0
691          ELSE
692             hicen(n) = vicen(n) / aicen(n)
693             hsnon(n) = vsnon(n) / aicen(n)
694             reduced_aicen(n) = c1 ! n=jpl
695             IF (n < jpl) reduced_aicen(n) = aicen(n) &
696                                  * (-0.024_wp*hicen(n) + 0.832_wp) 
697             asnon(n) = reduced_aicen(n)  ! effective snow fraction (empirical)
698             ! MV should check whether this makes sense to have the same effective snow fraction in here
699          END IF
700 
701 ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
702 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
703 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
704 ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
705 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
706 ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
707 ! alfan = 60% of the ice volume) in each category lies above the reference line,
708 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
709
710 ! MV: 
711 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE
712 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
713
714 ! Where does that choice come from
715 
716          alfan(n) = 0.6 * hicen(n)
717          betan(n) = 0.4 * hicen(n)
718       
719          cum_max_vol(n)     = z0
720          cum_max_vol_tmp(n) = z0
721     
722       END DO ! jpl
723 
724       cum_max_vol_tmp(0) = z0
725       drain = z0
726       dvolp = z0
727
728       !----------------------------------------------------------
729       ! x) Drain overflow water, update pond fraction and volume
730       !----------------------------------------------------------
731       
732       !--------------------------------------------------------------------------
733       ! the maximum amount of water that can be contained up to each ice category
734       !--------------------------------------------------------------------------
735
736       ! MV
737       ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
738       ! Then the excess volume cum_max_vol(jl) drains out of the system
739       ! It should be added to wfx_pnd
740       ! END MV
741     
742       DO n = 1, jpl-1 ! last category can not hold any volume
743 
744          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
745 
746             ! total volume in level including snow
747             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
748                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 
749 
750             ! subtract snow solid volumes from lower categories in current level
751             DO ns = 1, n 
752                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
753                   - rhosn/rhofw * &   ! free air fraction that can be filled by water
754                     asnon(ns)  * &    ! effective areal fraction of snow in that category
755                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
756                                   alfan(n)), z0) 
757             END DO
758
759          ELSE ! assume higher categories unoccupied
760             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
761          END IF
762          !IF (cum_max_vol_tmp(n) < z0) THEN
763          !   call abort_ice('negative melt pond volume')
764          !END IF
765       END DO
766       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
767       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
768     
769       !----------------------------------------------------------------
770       ! is there more meltwater than can be held in the floe?
771       !----------------------------------------------------------------
772       IF (zvolp >= cum_max_vol(jpl)) THEN
773          drain = zvolp - cum_max_vol(jpl) + epsi10
774          zvolp = zvolp - drain ! update meltwater volume available
775          dvolp = drain         ! this is the drained water
776          IF (zvolp < epsi10) THEN
777             dvolp = dvolp + zvolp
778             zvolp = z0
779          END IF
780       END IF
781     
782       ! height and area corresponding to the remaining volume
783 
784!      call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
785!           zvolp, cum_max_vol, hpond, m_index)
786     
787       DO n=1, m_index
788          zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update
789                                                   !  volume instead, no ?
790          zapondn(n) = reduced_aicen(n) 
791          ! in practise, pond fraction depends on the empirical snow fraction
792          ! so in turn on ice thickness
793       END DO
794     
795       !------------------------------------------------------------------------
796       ! Drainage through brine network (permeability)
797       !------------------------------------------------------------------------
798       !!! drainage due to ice permeability - Darcy's law
799     
800       ! sea water level
801       msno = z0
802       DO n=1,jpl
803         msno = msno + vsnon(n) * rhosn
804       END DO
805       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice
806       hsl_rel = floe_weight / rau0 &
807               - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
808     
809       deltah = hpond - hsl_rel
810       pressure_head = grav * rau0 * max(deltah, z0)
811 
812       ! drain IF ice is permeable   
813       permflag = 0
814       IF (pressure_head > z0) THEN
815       DO n = 1, jpl-1
816          IF (hicen(n) /= z0) THEN
817             perm = 0. ! MV ugly dummy patch
818             CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm)
819             IF (perm > z0) permflag = 1
820
821             drain = perm*zapondn(n)*pressure_head*rdt_ice / &
822                                      (viscosity*hicen(n))
823             dvolp = dvolp + min(drain, zvolp)
824             zvolp = max(zvolp - drain, z0)
825             IF (zvolp < epsi10) THEN
826                dvolp = dvolp + zvolp
827                zvolp = z0
828             END IF
829          END IF
830       END DO
831 
832       ! adjust melt pond DIMENSIONs
833       IF (permflag > 0) THEN
834          ! recompute pond depth   
835!         CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
836!                         zvolp, cum_max_vol, hpond, m_index)
837          DO n=1, m_index
838             zhpondn(n) = hpond - alfan(n) + alfan(1)
839             zapondn(n) = reduced_aicen(n) 
840          END DO
841       END IF
842       END IF ! pressure_head
843 
844       !-------------------------------
845       ! X) remove water from the snow
846       !-------------------------------
847       !------------------------------------------------------------------------
848       ! total melt pond volume in category DOes not include snow volume
849       ! snow in melt ponds is not melted
850       !------------------------------------------------------------------------
851 
852       ! Calculate pond volume for lower categories
853       DO n=1,m_index-1
854          zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow
855                   - (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n))
856       END DO
857 
858       ! Calculate pond volume for highest category = remaining pond volume
859
860       ! The following is completely unclear to Martin at least
861       ! Could we redefine properly and recode in a more readable way ?
862
863       ! m_index = last category with melt pond
864
865       IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water
866
867       IF (m_index > 1) THEN
868         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
869           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) !
870         ELSE
871           zvolpn(m_index) = z0
872           zhpondn(m_index) = z0
873           zapondn(m_index) = z0
874           ! If remaining pond volume is negative reduce pond volume of
875           ! lower category
876           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 
877             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
878                                + zvolp
879         END IF
880       END IF
881 
882       DO n=1,m_index
883          IF (zapondn(n) > epsi10) THEN
884              zhpondn(n) = zvolpn(n) / zapondn(n)
885          ELSE
886             dvolp = dvolp + zvolpn(n)
887             zhpondn(n) = z0
888             zvolpn(n) = z0
889             zapondn(n) = z0
890          end IF
891       END DO
892       DO n = m_index+1, jpl
893          zhpondn(n) = z0
894          zapondn(n) = z0
895          zvolpn (n) = z0
896       END DO
897 
898    END SUBROUTINE lim_mp_area
899
900!OLI_CODE   
901!OLI_CODE
902!OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
903!OLI_CODE                          zvolp, cum_max_vol,                &
904!OLI_CODE                          hpond, m_index)
905!OLI_CODE       !!-------------------------------------------------------------------
906!OLI_CODE       !!                ***  ROUTINE calc_hpond  ***
907!OLI_CODE       !!
908!OLI_CODE       !! ** Purpose :   Compute melt pond depth
909!OLI_CODE       !!-------------------------------------------------------------------
910!OLI_CODE     
911!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
912!OLI_CODE          aicen, &
913!OLI_CODE          asnon, &
914!OLI_CODE          hsnon, &
915!OLI_CODE          rhos,  &
916!OLI_CODE          alfan, &
917!OLI_CODE          cum_max_vol
918!OLI_CODE     
919!OLI_CODE       REAL (wp), INTENT(IN) :: &
920!OLI_CODE          zvolp
921!OLI_CODE     
922!OLI_CODE       REAL (wp), INTENT(OUT) :: &
923!OLI_CODE          hpond
924!OLI_CODE     
925!OLI_CODE       INTEGER, INTENT(OUT) :: &
926!OLI_CODE          m_index
927!OLI_CODE     
928!OLI_CODE       INTEGER :: n, ns
929!OLI_CODE     
930!OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: &
931!OLI_CODE          hitl, &
932!OLI_CODE          aicetl
933!OLI_CODE     
934!OLI_CODE       REAL (wp) :: &
935!OLI_CODE          rem_vol, &
936!OLI_CODE          area, &
937!OLI_CODE          vol, &
938!OLI_CODE          tmp, &
939!OLI_CODE          z0   = 0.0_wp,    &   
940!OLI_CODE          epsi10 = 1.0e-11_wp
941!OLI_CODE     
942!OLI_CODE       !----------------------------------------------------------------
943!OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained?
944!OLI_CODE       !----------------------------------------------------------------
945!OLI_CODE     
946!OLI_CODE       IF (zvolp < epsi10) THEN
947!OLI_CODE        hpond = z0
948!OLI_CODE        m_index = 0
949!OLI_CODE       ELSE
950!OLI_CODE       
951!OLI_CODE        !----------------------------------------------------------------
952!OLI_CODE        ! Calculate the category where water fills up to
953!OLI_CODE        !----------------------------------------------------------------
954!OLI_CODE       
955!OLI_CODE        !----------|
956!OLI_CODE        !          |
957!OLI_CODE        !          |
958!OLI_CODE        !          |----------|                                     -- --
959!OLI_CODE        !__________|__________|_________________________________________ ^
960!OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled
961!OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
962!OLI_CODE        !          |          |          |              |
963!OLI_CODE        !          |          |          |              |hpond
964!OLI_CODE        !          |          |          |----------|   |     |-------
965!OLI_CODE        !          |          |          |          |   |     |
966!OLI_CODE        !          |          |          |          |---v-----|
967!OLI_CODE        !          |          | m_index  |          |         |
968!OLI_CODE        !-------------------------------------------------------------
969!OLI_CODE       
970!OLI_CODE        m_index = 0  ! 1:m_index categories have water in them
971!OLI_CODE        DO n = 1, jpl
972!OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN
973!OLI_CODE              m_index = n
974!OLI_CODE              IF (n == 1) THEN
975!OLI_CODE                 rem_vol = zvolp
976!OLI_CODE              ELSE
977!OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1)
978!OLI_CODE              END IF
979!OLI_CODE              exit ! to break out of the loop
980!OLI_CODE           END IF
981!OLI_CODE        END DO
982!OLI_CODE        m_index = min(jpl-1, m_index)
983!OLI_CODE       
984!OLI_CODE        !----------------------------------------------------------------
985!OLI_CODE        ! semi-filled layer may have m_index different snow in it
986!OLI_CODE        !----------------------------------------------------------------
987!OLI_CODE       
988!OLI_CODE        !-----------------------------------------------------------  ^
989!OLI_CODE        !                                                             |  alfan(m_index+1)
990!OLI_CODE        !                                                             |
991!OLI_CODE        !hitl(3)-->                             |----------|          |
992!OLI_CODE        !hitl(2)-->                |------------| * * * * *|          |
993!OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
994!OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^
995!OLI_CODE        !                various snow from lower categories          |  |alfa(m_index)
996!OLI_CODE       
997!OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories
998!OLI_CODE        ! aicetl - area of each snow depth in this layer
999!OLI_CODE       
1000!OLI_CODE        hitl(:) = z0
1001!OLI_CODE        aicetl(:) = z0
1002!OLI_CODE        DO n = 1, m_index
1003!OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1004!OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0)
1005!OLI_CODE           aicetl(n) = asnon(n)
1006!OLI_CODE           
1007!OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1008!OLI_CODE        END DO
1009!OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1010!OLI_CODE        aicetl(m_index+1) = z0
1011!OLI_CODE       
1012!OLI_CODE        !----------------------------------------------------------------
1013!OLI_CODE        ! reorder array according to hitl
1014!OLI_CODE        ! snow heights not necessarily in height order
1015!OLI_CODE        !----------------------------------------------------------------
1016!OLI_CODE       
1017!OLI_CODE        DO ns = 1, m_index+1
1018!OLI_CODE           DO n = 0, m_index - ns + 1
1019!OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order
1020!OLI_CODE                 tmp = hitl(n)
1021!OLI_CODE                 hitl(n) = hitl(n+1)
1022!OLI_CODE                 hitl(n+1) = tmp
1023!OLI_CODE                 tmp = aicetl(n)
1024!OLI_CODE                 aicetl(n) = aicetl(n+1)
1025!OLI_CODE                 aicetl(n+1) = tmp
1026!OLI_CODE              END IF
1027!OLI_CODE           END DO
1028!OLI_CODE        END DO
1029!OLI_CODE       
1030!OLI_CODE        !----------------------------------------------------------------
1031!OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous
1032!OLI_CODE        !----------------------------------------------------------------
1033!OLI_CODE       
1034!OLI_CODE        !hitl(3)----------------------------------------------------------------
1035!OLI_CODE        !                                                       | * * * * * * * * 
1036!OLI_CODE        !                                                       |* * * * * * * * *
1037!OLI_CODE        !hitl(2)----------------------------------------------------------------
1038!OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * * 
1039!OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * *
1040!OLI_CODE        !hitl(1)----------------------------------------------------------------
1041!OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
1042!OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
1043!OLI_CODE        !hitl(0)----------------------------------------------------------------
1044!OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)           
1045!OLI_CODE       
1046!OLI_CODE        ! move up over layers incrementing volume
1047!OLI_CODE        DO n = 1, m_index+1
1048!OLI_CODE           
1049!OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer
1050!OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1051!OLI_CODE           
1052!OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
1053!OLI_CODE           
1054!OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
1055!OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
1056!OLI_CODE                           alfan(1)
1057!OLI_CODE              exit
1058!OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth
1059!OLI_CODE              rem_vol = rem_vol - vol
1060!OLI_CODE           END IF
1061!OLI_CODE           
1062!OLI_CODE        END DO
1063!OLI_CODE       
1064!OLI_CODE       END IF
1065!OLI_CODE     
1066!OLI_CODE    END SUBROUTINE calc_hpond
1067!OLI_CODE   
1068!OLI_CODE
1069    SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm)
1070       !!-------------------------------------------------------------------
1071       !!                ***  ROUTINE lim_mp_perm ***
1072       !!
1073       !! ** Purpose :   Determine the liquid fraction of brine in the ice
1074       !!                and its permeability
1075       !!-------------------------------------------------------------------
1076       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1077          ticen,  & ! energy of melting for each ice layer (J/m2)
1078          salin
1079 
1080       REAL (wp), INTENT(IN) :: &
1081          vicen     ! ice volume
1082     
1083       REAL (wp), INTENT(OUT) :: &
1084          perm      ! permeability
1085 
1086       REAL (wp) ::   &
1087          Sbr        ! brine salinity
1088 
1089       REAL (wp), DIMENSION(nlay_i) ::   &
1090          Tin, &    ! ice temperature
1091          phi       ! liquid fraction
1092 
1093       INTEGER :: k
1094     
1095       REAL (wp) :: &
1096          c2    = 2.0_wp
1097 
1098       !-----------------------------------------------------------------
1099       ! Compute ice temperatures from enthalpies using quadratic formula
1100       !-----------------------------------------------------------------
1101 
1102       DO k = 1,nlay_i
1103          Tin(k) = ticen(k) 
1104       END DO
1105 
1106       !-----------------------------------------------------------------
1107       ! brine salinity and liquid fraction
1108       !-----------------------------------------------------------------
1109 
1110       IF (maxval(Tin-rtt) <= -c2) THEN
1111 
1112          DO k = 1,nlay_i
1113             Sbr = - 1.2_wp                 &
1114                   -21.8_wp     * (Tin(k)-rtt)    &
1115                   - 0.919_wp   * (Tin(k)-rtt)**2 &
1116                   - 0.01878_wp * (Tin(k)-rtt)**3
1117             phi(k) = salin(k)/Sbr ! liquid fraction
1118          END DO ! k
1119       
1120       ELSE
1121 
1122          DO k = 1,nlay_i
1123             Sbr = -17.6_wp    * (Tin(k)-rtt)    &
1124                   - 0.389_wp  * (Tin(k)-rtt)**2 &
1125                   - 0.00362_wp* (Tin(k)-rtt)**3
1126             phi(k) = salin(k)/Sbr ! liquid fraction
1127          END DO
1128 
1129       END IF
1130     
1131       !-----------------------------------------------------------------
1132       ! permeability
1133       !-----------------------------------------------------------------
1134 
1135       perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking
1136                                            ! bastard of Golden)
1137     
1138    END SUBROUTINE lim_mp_perm
1139!OLI_CODE   
1140!OLI_CODE #else
1141!OLI_CODE    !!----------------------------------------------------------------------
1142!OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model
1143!OLI_CODE    !!----------------------------------------------------------------------
1144!OLI_CODE CONTAINS
1145!OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine
1146!OLI_CODE    END SUBROUTINE lim_mp_init   
1147!OLI_CODE    SUBROUTINE lim_mp               ! Empty routine
1148!OLI_CODE    END SUBROUTINE lim_mp
1149!OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine
1150!OLI_CODE    END SUBROUTINE compute_mp_topo       
1151!OLI_CODE    SUBROUTINE pond_area            ! Empty routine
1152!OLI_CODE    END SUBROUTINE pond_area   
1153!OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine
1154!OLI_CODE    END SUBROUTINE calc_hpond   
1155!OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine
1156!OLI_CODE    END SUBROUTINE permeability_phy   
1157!OLI_CODE #endif
1158!OLI_CODE    !!======================================================================
1159!OLI_CODE END MODULE limmp_topo
1160
1161
1162!OLI_CODE MODULE limmp_topo
1163!OLI_CODE    !!======================================================================
1164!OLI_CODE    !!                       ***  MODULE limmp_topo ***
1165!OLI_CODE    !! LIM-3 sea-ice :  computation of melt ponds' properties
1166!OLI_CODE    !!======================================================================
1167!OLI_CODE    !! History :  Original code by Daniela Flocco and Adrian Turner
1168!OLI_CODE    !!            ! 2012-09 (O. Lecomte) Adaptation for routine inclusion in
1169!OLI_CODE    !!                      NEMO-LIM3.1
1170!OLI_CODE    !!            ! 2016-11 (O. Lecomte, C. Rousset, M. Vancoppenolle)
1171!OLI_CODE    !!                      Adaptation for merge with NEMO-LIM3.6
1172!OLI_CODE    !!----------------------------------------------------------------------
1173!OLI_CODE #if defined key_lim3
1174!OLI_CODE    !!----------------------------------------------------------------------
1175!OLI_CODE    !!   'key_lim3'                                      LIM-3 sea-ice model
1176!OLI_CODE    !!----------------------------------------------------------------------
1177!OLI_CODE    !!   lim_mp_init     : melt pond properties initialization
1178!OLI_CODE    !!   lim_mp          : melt pond routine caller
1179!OLI_CODE    !!   compute_mp_topo : Actual melt pond routine
1180!OLI_CODE    !!----------------------------------------------------------------------
1181!OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice
1182!OLI_CODE    USE phycst
1183!OLI_CODE    USE dom_ice
1184!OLI_CODE    USE dom_oce
1185!OLI_CODE    USE sbc_oce
1186!OLI_CODE    USE sbc_ice
1187!OLI_CODE    USE par_ice
1188!OLI_CODE    USE par_oce
1189!OLI_CODE    USE ice
1190!OLI_CODE    USE thd_ice
1191!OLI_CODE    USE in_out_manager
1192!OLI_CODE    USE lbclnk
1193!OLI_CODE    USE lib_mpp
1194!OLI_CODE
1195!OLI_CODE    IMPLICIT NONE
1196!OLI_CODE    PRIVATE
1197!OLI_CODE
1198!OLI_CODE    PUBLIC   lim_mp_init
1199!OLI_CODE    PUBLIC   lim_mp
1200!OLI_CODE
1201!OLI_CODE CONTAINS
1202!OLI_CODE
1203!OLI_CODE    SUBROUTINE lim_mp_init
1204!OLI_CODE       !!-------------------------------------------------------------------
1205!OLI_CODE       !!                ***  ROUTINE lim_mp_init  ***
1206!OLI_CODE       !!
1207!OLI_CODE       !! ** Purpose :   Initialize melt ponds
1208!OLI_CODE       !!-------------------------------------------------------------------
1209!OLI_CODE       a_ip_frac(:,:,:)    = 0._wp
1210!OLI_CODE       a_ip(:,:,:)         = 0._wp
1211!OLI_CODE       h_ip(:,:,:)         = 0._wp
1212!OLI_CODE       v_ip(:,:,:)         = 0._wp
1213!OLI_CODE       h_il(:,:,:)         = 0._wp
1214!OLI_CODE       v_il(:,:,:)         = 0._wp
1215!OLI_CODE         
1216!OLI_CODE    END SUBROUTINE lim_mp_init
1217!OLI_CODE
1218!OLI_CODE
1219!OLI_CODE    SUBROUTINE lim_mp
1220!OLI_CODE       !!-------------------------------------------------------------------
1221!OLI_CODE       !!                ***  ROUTINE lim_mp  ***
1222!OLI_CODE       !!
1223!OLI_CODE       !! ** Purpose :   Compute surface heat flux and call main melt pond
1224!OLI_CODE       !!                routine
1225!OLI_CODE       !!-------------------------------------------------------------------
1226!OLI_CODE
1227!OLI_CODE       INTEGER  ::   ji, jj, jl   ! dummy loop indices
1228!OLI_CODE
1229!OLI_CODE       fsurf(:,:) = 0.e0
1230!OLI_CODE       DO jl = 1, jpl
1231!OLI_CODE          DO jj = 1, jpj
1232!OLI_CODE             DO ji = 1, jpi
1233!OLI_CODE                fsurf(ji,jj) = fsurf(ji,jj) + a_i(ji,jj,jl) * &
1234!OLI_CODE                      (qns_ice(ji,jj,jl) + (1.0 - izero(ji,jj,jl)) &
1235!OLI_CODE                        * qsr_ice(ji,jj,jl))
1236!OLI_CODE             END DO
1237!OLI_CODE          END DO
1238!OLI_CODE       END DO
1239!OLI_CODE
1240!OLI_CODE       CALL compute_mp_topo(at_i, a_i,                               &
1241!OLI_CODE                          vt_i, v_i, v_s, rhosn_glo, t_i, s_i, a_ip_frac,  &
1242!OLI_CODE                       h_ip, h_il, t_su, tatm_ice, diag_sur_me*rdt_ice, &
1243!OLI_CODE                          fsurf, fwoc)
1244!OLI_CODE
1245!OLI_CODE       at_ip(:,:) = 0.0
1246!OLI_CODE       vt_ip(:,:) = 0.0
1247!OLI_CODE       vt_il(:,:) = 0.0
1248!OLI_CODE       DO jl = 1, jpl
1249!OLI_CODE         DO jj = 1, jpj
1250!OLI_CODE            DO ji = 1, jpi
1251!OLI_CODE               a_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1252!OLI_CODE                                                   a_i(ji,jj,jl))
1253!OLI_CODE               v_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1254!OLI_CODE                                   a_i(ji,jj,jl) * h_ip(ji,jj,jl))
1255!OLI_CODE               v_il(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
1256!OLI_CODE                                   a_i(ji,jj,jl) * h_il(ji,jj,jl))
1257!OLI_CODE               at_ip(ji,jj)   = at_ip(ji,jj) + a_ip(ji,jj,jl)
1258!OLI_CODE               vt_ip(ji,jj)   = vt_ip(ji,jj) + v_ip(ji,jj,jl)
1259!OLI_CODE               vt_il(ji,jj)   = vt_il(ji,jj) + v_il(ji,jj,jl)
1260!OLI_CODE            END DO
1261!OLI_CODE         END DO
1262!OLI_CODE       END DO
1263!OLI_CODE
1264!OLI_CODE    END SUBROUTINE lim_mp
1265!OLI_CODE
1266!OLI_CODE
1267!OLI_CODE    SUBROUTINE compute_mp_topo(aice,      aicen,     &
1268!OLI_CODE                               vice,      vicen,     &
1269!OLI_CODE                               vsnon,     rhos,      &
1270!OLI_CODE                               ticen,     salin,     &
1271!OLI_CODE                               a_ip_frac, h_ip,      &
1272!OLI_CODE                               h_il,      Tsfc,      &
1273!OLI_CODE                               potT,      meltt,     &
1274!OLI_CODE                               fsurf,     fwoc)
1275!OLI_CODE       !!-------------------------------------------------------------------
1276!OLI_CODE       !!                ***  ROUTINE compute_mp_topo  ***
1277!OLI_CODE       !!
1278!OLI_CODE       !! ** Purpose :   Compute melt pond evolution based on the ice
1279!OLI_CODE       !!                topography as inferred from the ice thickness
1280!OLI_CODE       !!                distribution. 
1281!OLI_CODE       !!
1282!OLI_CODE       !! ** Method  :   This code is initially based on Flocco and Feltham
1283!OLI_CODE       !!                (2007) and Flocco et al. (2010). More to come...
1284!OLI_CODE       !!
1285!OLI_CODE       !! ** Tunable parameters :
1286!OLI_CODE       !!
1287!OLI_CODE       !! ** Note :
1288!OLI_CODE       !!
1289!OLI_CODE       !! ** References
1290!OLI_CODE       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond
1291!OLI_CODE       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:
1292!OLI_CODE       !!    10.1029/2006JC003836.
1293!OLI_CODE       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of
1294!OLI_CODE       !!    a physically based melt pond scheme into the sea ice component of a
1295!OLI_CODE       !!    climate model.  J. Geophys. Res. 115, C08012,
1296!OLI_CODE       !!    doi: 10.1029/2009JC005568.
1297!OLI_CODE       !!   
1298!OLI_CODE       !!-------------------------------------------------------------------
1299!OLI_CODE
1300!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
1301!OLI_CODE          INTENT(IN) :: &
1302!OLI_CODE          aice, &    ! total ice area fraction
1303!OLI_CODE          vice       ! total ice volume (m)
1304!OLI_CODE
1305!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
1306!OLI_CODE          INTENT(IN) :: &
1307!OLI_CODE          aicen, &   ! ice area fraction, per category
1308!OLI_CODE          vsnon, &   ! snow volume, per category (m)
1309!OLI_CODE          rhos,  &   ! equivalent snow density, per category (kg/m^3)
1310!OLI_CODE          vicen      ! ice volume, per category (m)
1311!OLI_CODE
1312!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
1313!OLI_CODE          INTENT(IN) :: &
1314!OLI_CODE          ticen, &   ! ice enthalpy, per category
1315!OLI_CODE          salin
1316!OLI_CODE
1317!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
1318!OLI_CODE          INTENT(INOUT) :: &
1319!OLI_CODE          a_ip_frac , &   ! pond area fraction of ice, per ice category
1320!OLI_CODE          h_ip , &   ! pond depth, per ice category (m)
1321!OLI_CODE          h_il       ! Refrozen ice lid thickness, per ice category (m)
1322!OLI_CODE
1323!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
1324!OLI_CODE          INTENT(IN) :: &
1325!OLI_CODE          potT,  &   ! air potential temperature
1326!OLI_CODE          meltt, &   ! total surface meltwater flux
1327!OLI_CODE          fsurf      ! thermodynamic heat flux at ice/snow surface (W/m^2)
1328!OLI_CODE
1329!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
1330!OLI_CODE          INTENT(INOUT) :: &
1331!OLI_CODE          fwoc      ! fresh water flux to the ocean (from draining and other pond volume adjustments)
1332!OLI_CODE                    ! (m)
1333!OLI_CODE
1334!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
1335!OLI_CODE          INTENT(IN) :: &
1336!OLI_CODE          Tsfc       ! snow/sea ice surface temperature
1337!OLI_CODE
1338!OLI_CODE       ! local variables
1339!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
1340!OLI_CODE          zTsfcn, & ! ice/snow surface temperature (C)
1341!OLI_CODE          zvolpn, & ! pond volume per unit area, per category (m)
1342!OLI_CODE          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
1343!OLI_CODE
1344!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
1345!OLI_CODE          zapondn,& ! pond area fraction, per category
1346!OLI_CODE          zhpondn   ! pond depth, per category (m)
1347!OLI_CODE
1348!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj) :: &
1349!OLI_CODE          zvolp       ! total volume of pond, per unit area of pond (m)
1350!OLI_CODE
1351!OLI_CODE       REAL (wp) :: &
1352!OLI_CODE          zhi,    & ! ice thickness (m)
1353!OLI_CODE          zdHui,  & ! change in thickness of ice lid (m)
1354!OLI_CODE          zomega, & ! conduction
1355!OLI_CODE          zdTice, & ! temperature difference across ice lid (C)
1356!OLI_CODE          zdvice, & ! change in ice volume (m)
1357!OLI_CODE          zTavg,  & ! mean surface temperature across categories (C)
1358!OLI_CODE          zTp,    & ! pond freezing temperature (C)
1359!OLI_CODE          zdvn      ! change in melt pond volume for fresh water budget
1360!OLI_CODE       INTEGER, DIMENSION (jpi*jpj) :: &
1361!OLI_CODE          indxi, indxj    ! compressed indices for cells with ice melting
1362!OLI_CODE
1363!OLI_CODE       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
1364!OLI_CODE
1365!OLI_CODE       INTEGER, DIMENSION (jpl) :: &
1366!OLI_CODE          kcells          ! cells where ice lid combines with vice
1367!OLI_CODE
1368!OLI_CODE       INTEGER, DIMENSION (jpi*jpj,jpl) :: &
1369!OLI_CODE          indxii, indxjj  ! i,j indices for kcells loop
1370!OLI_CODE
1371!OLI_CODE       REAL (wp), parameter :: &
1372!OLI_CODE          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m)
1373!OLI_CODE          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C)
1374!OLI_CODE          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
1375!OLI_CODE          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
1376!OLI_CODE          z0       = 0._wp,    &
1377!OLI_CODE          zTimelt   = 0._wp,    &
1378!OLI_CODE          z01      = 0.01_wp,  &
1379!OLI_CODE          z25      = 0.25_wp,  &
1380!OLI_CODE          z5       = 0.5_wp,   &
1381!OLI_CODE          epsi10     = 1.0e-11_wp
1382!OLI_CODE       !---------------------------------------------------------------
1383!OLI_CODE       ! initialize
1384!OLI_CODE       !---------------------------------------------------------------
1385!OLI_CODE
1386!OLI_CODE       DO j = 1, jpj
1387!OLI_CODE          DO i = 1, jpi
1388!OLI_CODE             zvolp(i,j) = z0
1389!OLI_CODE          END DO
1390!OLI_CODE       END DO
1391!OLI_CODE       DO n = 1, jpl
1392!OLI_CODE          DO j = 1, jpj
1393!OLI_CODE             DO i = 1, jpi
1394!OLI_CODE                ! load tracers
1395!OLI_CODE                zvolp(i,j) = zvolp(i,j) + h_ip(i,j,n) &
1396!OLI_CODE                                      * a_ip_frac(i,j,n) * aicen(i,j,n)
1397!OLI_CODE                zTsfcn(i,j,n) = Tsfc(i,j,n) - rtt ! convert in Celsius - Oli
1398!OLI_CODE                zvuin (i,j,n) = h_il(i,j,n) &
1399!OLI_CODE                             * a_ip_frac(i,j,n) * aicen(i,j,n)
1400!OLI_CODE
1401!OLI_CODE                zhpondn(i,j,n) = z0     ! pond depth, per category
1402!OLI_CODE                zapondn(i,j,n) = z0     ! pond area,  per category
1403!OLI_CODE             END DO
1404!OLI_CODE          END DO
1405!OLI_CODE          indxii(:,n) = 0
1406!OLI_CODE          indxjj(:,n) = 0
1407!OLI_CODE          kcells  (n) = 0
1408!OLI_CODE       END DO
1409!OLI_CODE
1410!OLI_CODE       ! The freezing temperature for meltponds is assumed slightly below 0C,
1411!OLI_CODE       ! as if meltponds had a little salt in them.  The salt budget is not
1412!OLI_CODE       ! altered for meltponds, but if it were then an actual pond freezing
1413!OLI_CODE       ! temperature could be computed.
1414!OLI_CODE
1415!OLI_CODE       zTp = zTimelt - zTd
1416!OLI_CODE
1417!OLI_CODE       !-----------------------------------------------------------------
1418!OLI_CODE       ! Identify grid cells with ponds
1419!OLI_CODE       !-----------------------------------------------------------------
1420!OLI_CODE
1421!OLI_CODE       icells = 0
1422!OLI_CODE       DO j = 1, jpj
1423!OLI_CODE       DO i = 1, jpi
1424!OLI_CODE          zhi = z0
1425!OLI_CODE          IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j)
1426!OLI_CODE          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
1427!OLI_CODE             zvolp(i,j) > zmin_volp*aice(i,j)) THEN
1428!OLI_CODE             icells = icells + 1
1429!OLI_CODE             indxi(icells) = i
1430!OLI_CODE             indxj(icells) = j
1431!OLI_CODE          ELSE  ! remove ponds on thin ice
1432!OLI_CODE             !fpond(i,j) = fpond(i,j) - zvolp(i,j)
1433!OLI_CODE             zvolpn(i,j,:) = z0
1434!OLI_CODE             zvuin (i,j,:) = z0
1435!OLI_CODE             zvolp (i,j) = z0
1436!OLI_CODE          END IF
1437!OLI_CODE       END DO                     ! i
1438!OLI_CODE       END DO                     ! j
1439!OLI_CODE
1440!OLI_CODE       DO ij = 1, icells
1441!OLI_CODE          i = indxi(ij)
1442!OLI_CODE          j = indxj(ij)
1443!OLI_CODE
1444!OLI_CODE          !--------------------------------------------------------------
1445!OLI_CODE          ! calculate pond area and depth
1446!OLI_CODE          !--------------------------------------------------------------
1447!OLI_CODE          CALL pond_area(aice(i,j),vice(i,j),rhos(i,j,:), &
1448!OLI_CODE                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
1449!OLI_CODE                    ticen(i,j,:,:), salin(i,j,:,:), &
1450!OLI_CODE                    zvolpn(i,j,:), zvolp(i,j), &
1451!OLI_CODE                    zapondn(i,j,:),zhpondn(i,j,:), zdvn)
1452!OLI_CODE
1453!OLI_CODE          fwoc(i,j) = fwoc(i,j) + zdvn ! -> Goes to fresh water budget
1454!OLI_CODE
1455!OLI_CODE          ! mean surface temperature
1456!OLI_CODE          zTavg = z0
1457!OLI_CODE          DO n = 1, jpl
1458!OLI_CODE             zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
1459!OLI_CODE          END DO
1460!OLI_CODE          zTavg = zTavg / aice(i,j)
1461!OLI_CODE
1462!OLI_CODE          DO n = 1, jpl-1
1463!OLI_CODE                       
1464!OLI_CODE             IF (zvuin(i,j,n) > epsi10) THEN
1465!OLI_CODE
1466!OLI_CODE          !----------------------------------------------------------------
1467!OLI_CODE          ! melting: floating upper ice layer melts in whole or part
1468!OLI_CODE          !----------------------------------------------------------------
1469!OLI_CODE !               IF (zTsfcn(i,j,n) > zTp) THEN
1470!OLI_CODE                IF (zTavg > zTp) THEN
1471!OLI_CODE
1472!OLI_CODE                   zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n))
1473!OLI_CODE                   IF (zdvice > epsi10) THEN
1474!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) - zdvice
1475!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice
1476!OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   + zdvice
1477!OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   + zdvice
1478!OLI_CODE                       
1479!OLI_CODE                      IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN
1480!OLI_CODE                         ! ice lid melted and category is pond covered
1481!OLI_CODE                         zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n)
1482!OLI_CODE                         !fwoc(i,j)   = fwoc(i,j)   + zvuin(i,j,n)
1483!OLI_CODE                         zvuin(i,j,n)  = z0
1484!OLI_CODE                      END IF
1485!OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
1486!OLI_CODE                   END IF
1487!OLI_CODE
1488!OLI_CODE          !----------------------------------------------------------------
1489!OLI_CODE          ! freezing: existing upper ice layer grows
1490!OLI_CODE          !----------------------------------------------------------------
1491!OLI_CODE                ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp
1492!OLI_CODE
1493!OLI_CODE                 ! dIFferential growth of base of surface floating ice layer
1494!OLI_CODE                   zdTice = max(-zTavg, z0) ! > 0
1495!OLI_CODE                   zomega = rcdic*zdTice * zr1_rlfus
1496!OLI_CODE                   zdHui = sqrt(zomega*rdt_ice + z25*(zvuin(i,j,n)/  &
1497!OLI_CODE                         aicen(i,j,n))**2)- z5 * zvuin(i,j,n)/aicen(i,j,n)
1498!OLI_CODE
1499!OLI_CODE                   zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))   
1500!OLI_CODE                   IF (zdvice > epsi10) THEN
1501!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) + zdvice
1502!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
1503!OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   - zdvice
1504!OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   - zdvice
1505!OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
1506!OLI_CODE                   END IF
1507!OLI_CODE
1508!OLI_CODE                END IF ! zTavg
1509!OLI_CODE
1510!OLI_CODE          !----------------------------------------------------------------
1511!OLI_CODE          ! freezing: upper ice layer begins to form
1512!OLI_CODE          ! note: albedo does not change
1513!OLI_CODE          !----------------------------------------------------------------
1514!OLI_CODE             ELSE ! zvuin < epsi10
1515!OLI_CODE                     
1516!OLI_CODE                ! thickness of newly formed ice
1517!OLI_CODE                ! the surface temperature of a meltpond is the same as that
1518!OLI_CODE                ! of the ice underneath (0C), and the thermodynamic surface
1519!OLI_CODE                ! flux is the same
1520!OLI_CODE                zdHui = max(-fsurf(i,j)*rdt_ice*zr1_rlfus, z0)
1521!OLI_CODE                zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n)) 
1522!OLI_CODE                IF (zdvice > epsi10) THEN
1523!OLI_CODE                   zvuin (i,j,n) = zdvice
1524!OLI_CODE                   zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
1525!OLI_CODE                   zvolp (i,j)   = zvolp (i,j)   - zdvice
1526!OLI_CODE                   !fwoc(i,j)   = fwoc(i,j)   - zdvice
1527!OLI_CODE                   zhpondn(i,j,n)= zvolpn(i,j,n) / zapondn(i,j,n)
1528!OLI_CODE                END IF
1529!OLI_CODE                     
1530!OLI_CODE             END IF  ! zvuin
1531!OLI_CODE
1532!OLI_CODE          END DO ! jpl
1533!OLI_CODE
1534!OLI_CODE       END DO ! ij
1535!OLI_CODE
1536!OLI_CODE       !---------------------------------------------------------------
1537!OLI_CODE       ! remove ice lid if there is no liquid pond
1538!OLI_CODE       ! zvuin may be nonzero on category jpl due to dynamics
1539!OLI_CODE       !---------------------------------------------------------------
1540!OLI_CODE
1541!OLI_CODE       DO j = 1, jpj
1542!OLI_CODE       DO i = 1, jpi
1543!OLI_CODE          DO n = 1, jpl
1544!OLI_CODE             IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny &
1545!OLI_CODE                                     .and. zvuin (i,j,n) > epsi10) THEN
1546!OLI_CODE                kcells(n) = kcells(n) + 1
1547!OLI_CODE                indxij    = kcells(n)
1548!OLI_CODE                indxii(indxij,n) = i
1549!OLI_CODE                indxjj(indxij,n) = j
1550!OLI_CODE             END IF
1551!OLI_CODE          END DO
1552!OLI_CODE       END DO                     ! i
1553!OLI_CODE       END DO                     ! j
1554!OLI_CODE
1555!OLI_CODE       DO n = 1, jpl
1556!OLI_CODE
1557!OLI_CODE          IF (kcells(n) > 0) THEN
1558!OLI_CODE          DO ij = 1, kcells(n)
1559!OLI_CODE             i = indxii(ij,n)
1560!OLI_CODE             j = indxjj(ij,n)
1561!OLI_CODE             fwoc(i,j) = fwoc(i,j) + rhoic/rauw * zvuin(i,j,n) ! Completely refrozen lid goes into ocean (to be changed)
1562!OLI_CODE             zvuin(i,j,n) = z0
1563!OLI_CODE          END DO    ! ij
1564!OLI_CODE          END IF
1565!OLI_CODE
1566!OLI_CODE          ! reload tracers
1567!OLI_CODE          DO j = 1, jpj
1568!OLI_CODE             DO i = 1, jpi
1569!OLI_CODE                IF (zapondn(i,j,n) > epsi10) THEN
1570!OLI_CODE                   h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n)
1571!OLI_CODE                ELSE
1572!OLI_CODE                   zvuin(i,j,n) = z0
1573!OLI_CODE                   h_il(i,j,n) = z0
1574!OLI_CODE                END IF
1575!OLI_CODE                IF (aicen(i,j,n) > epsi10) THEN
1576!OLI_CODE                   a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * &
1577!OLI_CODE                         (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n))))
1578!OLI_CODE                   h_ip(i,j,n) = zhpondn(i,j,n)
1579!OLI_CODE                ELSE
1580!OLI_CODE                   a_ip_frac(i,j,n) = z0
1581!OLI_CODE                   h_ip(i,j,n) = z0
1582!OLI_CODE                   h_il(i,j,n) = z0
1583!OLI_CODE                END IF
1584!OLI_CODE             END DO ! i
1585!OLI_CODE          END DO    ! j
1586!OLI_CODE
1587!OLI_CODE       END DO       ! n
1588!OLI_CODE
1589!OLI_CODE    END SUBROUTINE compute_mp_topo
1590!OLI_CODE
1591!OLI_CODE
1592!OLI_CODE    SUBROUTINE pond_area(aice,vice,rhos,             &
1593!OLI_CODE                         aicen, vicen, vsnon, ticen, &
1594!OLI_CODE                         salin, zvolpn, zvolp,         &
1595!OLI_CODE                         zapondn,zhpondn,dvolp)
1596!OLI_CODE       !!-------------------------------------------------------------------
1597!OLI_CODE       !!                ***  ROUTINE pond_area  ***
1598!OLI_CODE       !!
1599!OLI_CODE       !! ** Purpose :   Compute melt pond area, depth and melting rates
1600!OLI_CODE       !!------------------------------------------------------------------
1601!OLI_CODE       REAL (wp), INTENT(IN) :: &
1602!OLI_CODE          aice,vice
1603!OLI_CODE
1604!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1605!OLI_CODE          aicen, vicen, vsnon, rhos
1606!OLI_CODE
1607!OLI_CODE       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
1608!OLI_CODE          ticen, salin
1609!OLI_CODE
1610!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
1611!OLI_CODE          zvolpn
1612!OLI_CODE
1613!OLI_CODE       REAL (wp), INTENT(INOUT) :: &
1614!OLI_CODE          zvolp, dvolp
1615!OLI_CODE
1616!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
1617!OLI_CODE          zapondn, zhpondn
1618!OLI_CODE
1619!OLI_CODE       INTEGER :: &
1620!OLI_CODE          n, ns,   &
1621!OLI_CODE          m_index, &
1622!OLI_CODE          permflag
1623!OLI_CODE
1624!OLI_CODE       REAL (wp), DIMENSION(jpl) :: &
1625!OLI_CODE          hicen, &
1626!OLI_CODE          hsnon, &
1627!OLI_CODE          asnon, &
1628!OLI_CODE          alfan, &
1629!OLI_CODE          betan, &
1630!OLI_CODE          cum_max_vol, &
1631!OLI_CODE          reduced_aicen       
1632!OLI_CODE
1633!OLI_CODE       REAL (wp), DIMENSION(0:jpl) :: &
1634!OLI_CODE          cum_max_vol_tmp
1635!OLI_CODE
1636!OLI_CODE       REAL (wp) :: &
1637!OLI_CODE          hpond, &
1638!OLI_CODE          drain, &
1639!OLI_CODE          floe_weight, &
1640!OLI_CODE          pressure_head, &
1641!OLI_CODE          hsl_rel, &
1642!OLI_CODE          deltah, &
1643!OLI_CODE          perm, &
1644!OLI_CODE          apond, &
1645!OLI_CODE          msno
1646!OLI_CODE
1647!OLI_CODE       REAL (wp), parameter :: &
1648!OLI_CODE          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s
1649!OLI_CODE          z0        = 0.0_wp    , &
1650!OLI_CODE          c1        = 1.0_wp    , &
1651!OLI_CODE          p4        = 0.4_wp    , &
1652!OLI_CODE          p6        = 0.6_wp    , &
1653!OLI_CODE          epsi10      = 1.0e-11_wp
1654!OLI_CODE         
1655!OLI_CODE      !-----------|
1656!OLI_CODE      !           |
1657!OLI_CODE      !           |-----------|
1658!OLI_CODE      !___________|___________|______________________________________sea-level
1659!OLI_CODE      !           |           |
1660!OLI_CODE      !           |           |---^--------|
1661!OLI_CODE      !           |           |   |        |
1662!OLI_CODE      !           |           |   |        |-----------|              |-------
1663!OLI_CODE      !           |           |   |alfan(n)|           |              |
1664!OLI_CODE      !           |           |   |        |           |--------------|
1665!OLI_CODE      !           |           |   |        |           |              |
1666!OLI_CODE      !---------------------------v-------------------------------------------
1667!OLI_CODE      !           |           |   ^        |           |              |
1668!OLI_CODE      !           |           |   |        |           |--------------|
1669!OLI_CODE      !           |           |   |betan(n)|           |              |
1670!OLI_CODE      !           |           |   |        |-----------|              |-------
1671!OLI_CODE      !           |           |   |        |
1672!OLI_CODE      !           |           |---v------- |
1673!OLI_CODE      !           |           |
1674!OLI_CODE      !           |-----------|
1675!OLI_CODE      !           |
1676!OLI_CODE      !-----------|
1677!OLI_CODE     
1678!OLI_CODE       !-------------------------------------------------------------------
1679!OLI_CODE       ! initialize
1680!OLI_CODE       !-------------------------------------------------------------------
1681!OLI_CODE
1682!OLI_CODE       DO n = 1, jpl
1683!OLI_CODE
1684!OLI_CODE          zapondn(n) = z0
1685!OLI_CODE          zhpondn(n) = z0
1686!OLI_CODE
1687!OLI_CODE          IF (aicen(n) < epsi10)  THEN
1688!OLI_CODE             hicen(n) =  z0
1689!OLI_CODE             hsnon(n) = z0
1690!OLI_CODE             reduced_aicen(n) = z0
1691!OLI_CODE          ELSE
1692!OLI_CODE             hicen(n) = vicen(n) / aicen(n)
1693!OLI_CODE             hsnon(n) = vsnon(n) / aicen(n)
1694!OLI_CODE             reduced_aicen(n) = c1 ! n=jpl
1695!OLI_CODE             IF (n < jpl) reduced_aicen(n) = aicen(n) &
1696!OLI_CODE                                  * (-0.024_wp*hicen(n) + 0.832_wp)
1697!OLI_CODE             asnon(n) = reduced_aicen(n)
1698!OLI_CODE          END IF
1699!OLI_CODE
1700!OLI_CODE ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
1701!OLI_CODE ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
1702!OLI_CODE ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
1703!OLI_CODE ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
1704!OLI_CODE ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
1705!OLI_CODE ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
1706!OLI_CODE ! alfan = 60% of the ice volume) in each category lies above the reference line,
1707!OLI_CODE ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
1708!OLI_CODE
1709!OLI_CODE          alfan(n) = p6 * hicen(n)
1710!OLI_CODE          betan(n) = p4 * hicen(n)
1711!OLI_CODE       
1712!OLI_CODE          cum_max_vol(n)     = z0
1713!OLI_CODE          cum_max_vol_tmp(n) = z0
1714!OLI_CODE     
1715!OLI_CODE       END DO ! jpl
1716!OLI_CODE
1717!OLI_CODE       cum_max_vol_tmp(0) = z0
1718!OLI_CODE       drain = z0
1719!OLI_CODE       dvolp = z0
1720!OLI_CODE     
1721!OLI_CODE       !--------------------------------------------------------------------------
1722!OLI_CODE       ! the maximum amount of water that can be contained up to each ice category
1723!OLI_CODE       !--------------------------------------------------------------------------
1724!OLI_CODE     
1725!OLI_CODE       DO n = 1, jpl-1 ! last category can not hold any volume
1726!OLI_CODE
1727!OLI_CODE          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
1728!OLI_CODE
1729!OLI_CODE             ! total volume in level including snow
1730!OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
1731!OLI_CODE                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
1732!OLI_CODE
1733!OLI_CODE
1734!OLI_CODE             ! subtract snow solid volumes from lower categories in current level
1735!OLI_CODE             DO ns = 1, n
1736!OLI_CODE                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
1737!OLI_CODE                   - rhos(ns)/rauw  * &    ! fraction of snow that is occupied by solid ??rauw
1738!OLI_CODE                     asnon(ns)  * &    ! area of snow from that category
1739!OLI_CODE                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
1740!OLI_CODE                                   alfan(n)), z0) 
1741!OLI_CODE                                       ! thickness of snow from ns layer in n layer
1742!OLI_CODE             END DO
1743!OLI_CODE
1744!OLI_CODE          ELSE ! assume higher categories unoccupied
1745!OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
1746!OLI_CODE          END IF
1747!OLI_CODE          !IF (cum_max_vol_tmp(n) < z0) THEN
1748!OLI_CODE          !   call abort_ice('negative melt pond volume')
1749!OLI_CODE          !END IF
1750!OLI_CODE       END DO
1751!OLI_CODE       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
1752!OLI_CODE       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
1753!OLI_CODE     
1754!OLI_CODE       !----------------------------------------------------------------
1755!OLI_CODE       ! is there more meltwater than can be held in the floe?
1756!OLI_CODE       !----------------------------------------------------------------
1757!OLI_CODE       IF (zvolp >= cum_max_vol(jpl)) THEN
1758!OLI_CODE          drain = zvolp - cum_max_vol(jpl) + epsi10
1759!OLI_CODE          zvolp = zvolp - drain
1760!OLI_CODE          dvolp = drain
1761!OLI_CODE          IF (zvolp < epsi10) THEN
1762!OLI_CODE             dvolp = dvolp + zvolp
1763!OLI_CODE             zvolp = z0
1764!OLI_CODE          END IF
1765!OLI_CODE       END IF
1766!OLI_CODE     
1767!OLI_CODE       ! height and area corresponding to the remaining volume
1768!OLI_CODE
1769!OLI_CODE       call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
1770!OLI_CODE            zvolp, cum_max_vol, hpond, m_index)
1771!OLI_CODE     
1772!OLI_CODE       DO n=1, m_index
1773!OLI_CODE          zhpondn(n) = hpond - alfan(n) + alfan(1)
1774!OLI_CODE          zapondn(n) = reduced_aicen(n)
1775!OLI_CODE       END DO
1776!OLI_CODE       apond = sum(zapondn(1:m_index))
1777!OLI_CODE     
1778!OLI_CODE       !------------------------------------------------------------------------
1779!OLI_CODE       ! drainage due to ice permeability - Darcy's law
1780!OLI_CODE       !------------------------------------------------------------------------
1781!OLI_CODE     
1782!OLI_CODE       ! sea water level
1783!OLI_CODE       msno = z0
1784!OLI_CODE       DO n=1,jpl
1785!OLI_CODE         msno = msno + vsnon(n) * rhos(n)
1786!OLI_CODE       END DO
1787!OLI_CODE       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice
1788!OLI_CODE       hsl_rel = floe_weight / rau0 &
1789!OLI_CODE               - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
1790!OLI_CODE     
1791!OLI_CODE       deltah = hpond - hsl_rel
1792!OLI_CODE       pressure_head = grav * rau0 * max(deltah, z0)
1793!OLI_CODE
1794!OLI_CODE       ! drain IF ice is permeable   
1795!OLI_CODE       permflag = 0
1796!OLI_CODE       IF (pressure_head > z0) THEN
1797!OLI_CODE       DO n = 1, jpl-1
1798!OLI_CODE          IF (hicen(n) /= z0) THEN
1799!OLI_CODE             CALL permeability_phi(ticen(:,n), salin(:,n), vicen(n), perm)
1800!OLI_CODE             IF (perm > z0) permflag = 1
1801!OLI_CODE             drain = perm*zapondn(n)*pressure_head*rdt_ice / &
1802!OLI_CODE                                      (viscosity*hicen(n))
1803!OLI_CODE             dvolp = dvolp + min(drain, zvolp)
1804!OLI_CODE             zvolp = max(zvolp - drain, z0)
1805!OLI_CODE             IF (zvolp < epsi10) THEN
1806!OLI_CODE                dvolp = dvolp + zvolp
1807!OLI_CODE                zvolp = z0
1808!OLI_CODE             END IF
1809!OLI_CODE          END IF
1810!OLI_CODE       END DO
1811!OLI_CODE 
1812!OLI_CODE       ! adjust melt pond DIMENSIONs
1813!OLI_CODE       IF (permflag > 0) THEN
1814!OLI_CODE          ! recompute pond depth   
1815!OLI_CODE          CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
1816!OLI_CODE                          zvolp, cum_max_vol, hpond, m_index)
1817!OLI_CODE          DO n=1, m_index
1818!OLI_CODE             zhpondn(n) = hpond - alfan(n) + alfan(1)
1819!OLI_CODE             zapondn(n) = reduced_aicen(n)
1820!OLI_CODE          END DO
1821!OLI_CODE          apond = sum(zapondn(1:m_index))
1822!OLI_CODE       END IF
1823!OLI_CODE       END IF ! pressure_head
1824!OLI_CODE
1825!OLI_CODE       !------------------------------------------------------------------------
1826!OLI_CODE       ! total melt pond volume in category DOes not include snow volume
1827!OLI_CODE       ! snow in melt ponds is not melted
1828!OLI_CODE       !------------------------------------------------------------------------
1829!OLI_CODE
1830!OLI_CODE       ! Calculate pond volume for lower categories
1831!OLI_CODE       DO n=1,m_index-1
1832!OLI_CODE          zvolpn(n) = zapondn(n) * zhpondn(n) &
1833!OLI_CODE                   - (rhos(n)/rauw) * asnon(n) * min(hsnon(n), zhpondn(n))!
1834!OLI_CODE       END DO
1835!OLI_CODE
1836!OLI_CODE       ! Calculate pond volume for highest category = remaining pond volume
1837!OLI_CODE       IF (m_index == 1) zvolpn(m_index) = zvolp
1838!OLI_CODE       IF (m_index > 1) THEN
1839!OLI_CODE         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
1840!OLI_CODE           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1))
1841!OLI_CODE         ELSE
1842!OLI_CODE           zvolpn(m_index) = z0
1843!OLI_CODE           zhpondn(m_index) = z0
1844!OLI_CODE           zapondn(m_index) = z0
1845!OLI_CODE           ! If remaining pond volume is negative reduce pond volume of
1846!OLI_CODE           ! lower category
1847!OLI_CODE           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &
1848!OLI_CODE             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
1849!OLI_CODE                                + zvolp
1850!OLI_CODE         END IF
1851!OLI_CODE       END IF
1852!OLI_CODE
1853!OLI_CODE       DO n=1,m_index
1854!OLI_CODE          IF (zapondn(n) > epsi10) THEN
1855!OLI_CODE              zhpondn(n) = zvolpn(n) / zapondn(n)
1856!OLI_CODE          ELSE
1857!OLI_CODE             dvolp = dvolp + zvolpn(n)
1858!OLI_CODE             zhpondn(n) = z0
1859!OLI_CODE             zvolpn(n) = z0
1860!OLI_CODE             zapondn(n) = z0
1861!OLI_CODE          end IF
1862!OLI_CODE       END DO
1863!OLI_CODE       DO n = m_index+1, jpl
1864!OLI_CODE          zhpondn(n) = z0
1865!OLI_CODE          zapondn(n) = z0
1866!OLI_CODE          zvolpn (n) = z0
1867!OLI_CODE       END DO
1868!OLI_CODE
1869!OLI_CODE    END SUBROUTINE pond_area
1870!OLI_CODE   
1871!OLI_CODE
1872!OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
1873!OLI_CODE                          zvolp, cum_max_vol,                &
1874!OLI_CODE                          hpond, m_index)
1875!OLI_CODE       !!-------------------------------------------------------------------
1876!OLI_CODE       !!                ***  ROUTINE calc_hpond  ***
1877!OLI_CODE       !!
1878!OLI_CODE       !! ** Purpose :   Compute melt pond depth
1879!OLI_CODE       !!-------------------------------------------------------------------
1880!OLI_CODE     
1881!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1882!OLI_CODE          aicen, &
1883!OLI_CODE          asnon, &
1884!OLI_CODE          hsnon, &
1885!OLI_CODE          rhos,  &
1886!OLI_CODE          alfan, &
1887!OLI_CODE          cum_max_vol
1888!OLI_CODE     
1889!OLI_CODE       REAL (wp), INTENT(IN) :: &
1890!OLI_CODE          zvolp
1891!OLI_CODE     
1892!OLI_CODE       REAL (wp), INTENT(OUT) :: &
1893!OLI_CODE          hpond
1894!OLI_CODE     
1895!OLI_CODE       INTEGER, INTENT(OUT) :: &
1896!OLI_CODE          m_index
1897!OLI_CODE     
1898!OLI_CODE       INTEGER :: n, ns
1899!OLI_CODE     
1900!OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: &
1901!OLI_CODE          hitl, &
1902!OLI_CODE          aicetl
1903!OLI_CODE     
1904!OLI_CODE       REAL (wp) :: &
1905!OLI_CODE          rem_vol, &
1906!OLI_CODE          area, &
1907!OLI_CODE          vol, &
1908!OLI_CODE          tmp, &
1909!OLI_CODE          z0   = 0.0_wp,    &   
1910!OLI_CODE          epsi10 = 1.0e-11_wp
1911!OLI_CODE     
1912!OLI_CODE       !----------------------------------------------------------------
1913!OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained?
1914!OLI_CODE       !----------------------------------------------------------------
1915!OLI_CODE     
1916!OLI_CODE       IF (zvolp < epsi10) THEN
1917!OLI_CODE        hpond = z0
1918!OLI_CODE        m_index = 0
1919!OLI_CODE       ELSE
1920!OLI_CODE       
1921!OLI_CODE        !----------------------------------------------------------------
1922!OLI_CODE        ! Calculate the category where water fills up to
1923!OLI_CODE        !----------------------------------------------------------------
1924!OLI_CODE       
1925!OLI_CODE        !----------|
1926!OLI_CODE        !          |
1927!OLI_CODE        !          |
1928!OLI_CODE        !          |----------|                                     -- --
1929!OLI_CODE        !__________|__________|_________________________________________ ^
1930!OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled
1931!OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
1932!OLI_CODE        !          |          |          |              |
1933!OLI_CODE        !          |          |          |              |hpond
1934!OLI_CODE        !          |          |          |----------|   |     |-------
1935!OLI_CODE        !          |          |          |          |   |     |
1936!OLI_CODE        !          |          |          |          |---v-----|
1937!OLI_CODE        !          |          | m_index  |          |         |
1938!OLI_CODE        !-------------------------------------------------------------
1939!OLI_CODE       
1940!OLI_CODE        m_index = 0  ! 1:m_index categories have water in them
1941!OLI_CODE        DO n = 1, jpl
1942!OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN
1943!OLI_CODE              m_index = n
1944!OLI_CODE              IF (n == 1) THEN
1945!OLI_CODE                 rem_vol = zvolp
1946!OLI_CODE              ELSE
1947!OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1)
1948!OLI_CODE              END IF
1949!OLI_CODE              exit ! to break out of the loop
1950!OLI_CODE           END IF
1951!OLI_CODE        END DO
1952!OLI_CODE        m_index = min(jpl-1, m_index)
1953!OLI_CODE       
1954!OLI_CODE        !----------------------------------------------------------------
1955!OLI_CODE        ! semi-filled layer may have m_index different snow in it
1956!OLI_CODE        !----------------------------------------------------------------
1957!OLI_CODE       
1958!OLI_CODE        !-----------------------------------------------------------  ^
1959!OLI_CODE        !                                                             |  alfan(m_index+1)
1960!OLI_CODE        !                                                             |
1961!OLI_CODE        !hitl(3)-->                             |----------|          |
1962!OLI_CODE        !hitl(2)-->                |------------| * * * * *|          |
1963!OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
1964!OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^
1965!OLI_CODE        !                various snow from lower categories          |  |alfa(m_index)
1966!OLI_CODE       
1967!OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories
1968!OLI_CODE        ! aicetl - area of each snow depth in this layer
1969!OLI_CODE       
1970!OLI_CODE        hitl(:) = z0
1971!OLI_CODE        aicetl(:) = z0
1972!OLI_CODE        DO n = 1, m_index
1973!OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1974!OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0)
1975!OLI_CODE           aicetl(n) = asnon(n)
1976!OLI_CODE           
1977!OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1978!OLI_CODE        END DO
1979!OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1980!OLI_CODE        aicetl(m_index+1) = z0
1981!OLI_CODE       
1982!OLI_CODE        !----------------------------------------------------------------
1983!OLI_CODE        ! reorder array according to hitl
1984!OLI_CODE        ! snow heights not necessarily in height order
1985!OLI_CODE        !----------------------------------------------------------------
1986!OLI_CODE       
1987!OLI_CODE        DO ns = 1, m_index+1
1988!OLI_CODE           DO n = 0, m_index - ns + 1
1989!OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order
1990!OLI_CODE                 tmp = hitl(n)
1991!OLI_CODE                 hitl(n) = hitl(n+1)
1992!OLI_CODE                 hitl(n+1) = tmp
1993!OLI_CODE                 tmp = aicetl(n)
1994!OLI_CODE                 aicetl(n) = aicetl(n+1)
1995!OLI_CODE                 aicetl(n+1) = tmp
1996!OLI_CODE              END IF
1997!OLI_CODE           END DO
1998!OLI_CODE        END DO
1999!OLI_CODE       
2000!OLI_CODE        !----------------------------------------------------------------
2001!OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous
2002!OLI_CODE        !----------------------------------------------------------------
2003!OLI_CODE       
2004!OLI_CODE        !hitl(3)----------------------------------------------------------------
2005!OLI_CODE        !                                                       | * * * * * * * * 
2006!OLI_CODE        !                                                       |* * * * * * * * *
2007!OLI_CODE        !hitl(2)----------------------------------------------------------------
2008!OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * * 
2009!OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * *
2010!OLI_CODE        !hitl(1)----------------------------------------------------------------
2011!OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
2012!OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
2013!OLI_CODE        !hitl(0)----------------------------------------------------------------
2014!OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)           
2015!OLI_CODE       
2016!OLI_CODE        ! move up over layers incrementing volume
2017!OLI_CODE        DO n = 1, m_index+1
2018!OLI_CODE           
2019!OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer
2020!OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
2021!OLI_CODE           
2022!OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
2023!OLI_CODE           
2024!OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
2025!OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
2026!OLI_CODE                           alfan(1)
2027!OLI_CODE              exit
2028!OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth
2029!OLI_CODE              rem_vol = rem_vol - vol
2030!OLI_CODE           END IF
2031!OLI_CODE           
2032!OLI_CODE        END DO
2033!OLI_CODE       
2034!OLI_CODE       END IF
2035!OLI_CODE     
2036!OLI_CODE    END SUBROUTINE calc_hpond
2037!OLI_CODE   
2038!OLI_CODE
2039!OLI_CODE    SUBROUTINE permeability_phi(ticen, salin, vicen, perm)
2040!OLI_CODE       !!-------------------------------------------------------------------
2041!OLI_CODE       !!                ***  ROUTINE permeability_phi  ***
2042!OLI_CODE       !!
2043!OLI_CODE       !! ** Purpose :   Determine the liquid fraction of brine in the ice
2044!OLI_CODE       !!                and its permeability
2045!OLI_CODE       !!-------------------------------------------------------------------
2046!OLI_CODE       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
2047!OLI_CODE          ticen,  & ! energy of melting for each ice layer (J/m2)
2048!OLI_CODE          salin
2049!OLI_CODE
2050!OLI_CODE       REAL (wp), INTENT(IN) :: &
2051!OLI_CODE          vicen     ! ice volume
2052!OLI_CODE     
2053!OLI_CODE       REAL (wp), INTENT(OUT) :: &
2054!OLI_CODE          perm      ! permeability
2055!OLI_CODE
2056!OLI_CODE       REAL (wp) ::   &
2057!OLI_CODE          Sbr        ! brine salinity
2058!OLI_CODE
2059!OLI_CODE       REAL (wp), DIMENSION(nlay_i) ::   &
2060!OLI_CODE          Tin, &    ! ice temperature
2061!OLI_CODE          phi       ! liquid fraction
2062!OLI_CODE
2063!OLI_CODE       INTEGER :: k
2064!OLI_CODE     
2065!OLI_CODE       REAL (wp) :: &
2066!OLI_CODE          c2    = 2.0_wp
2067!OLI_CODE 
2068!OLI_CODE       !-----------------------------------------------------------------
2069!OLI_CODE       ! Compute ice temperatures from enthalpies using quadratic formula
2070!OLI_CODE       !-----------------------------------------------------------------
2071!OLI_CODE
2072!OLI_CODE       DO k = 1,nlay_i
2073!OLI_CODE          Tin(k) = ticen(k)
2074!OLI_CODE       END DO
2075!OLI_CODE
2076!OLI_CODE       !-----------------------------------------------------------------
2077!OLI_CODE       ! brine salinity and liquid fraction
2078!OLI_CODE       !-----------------------------------------------------------------
2079!OLI_CODE
2080!OLI_CODE       IF (maxval(Tin-rtt) <= -c2) THEN
2081!OLI_CODE
2082!OLI_CODE          DO k = 1,nlay_i
2083!OLI_CODE             Sbr = - 1.2_wp                 &
2084!OLI_CODE                   -21.8_wp     * (Tin(k)-rtt)    &
2085!OLI_CODE                   - 0.919_wp   * (Tin(k)-rtt)**2 &
2086!OLI_CODE                   - 0.01878_wp * (Tin(k)-rtt)**3
2087!OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction
2088!OLI_CODE          END DO ! k
2089!OLI_CODE       
2090!OLI_CODE       ELSE
2091!OLI_CODE
2092!OLI_CODE          DO k = 1,nlay_i
2093!OLI_CODE             Sbr = -17.6_wp    * (Tin(k)-rtt)    &
2094!OLI_CODE                   - 0.389_wp  * (Tin(k)-rtt)**2 &
2095!OLI_CODE                   - 0.00362_wp* (Tin(k)-rtt)**3
2096!OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction
2097!OLI_CODE          END DO
2098!OLI_CODE
2099!OLI_CODE       END IF
2100!OLI_CODE     
2101!OLI_CODE       !-----------------------------------------------------------------
2102!OLI_CODE       ! permeability
2103!OLI_CODE       !-----------------------------------------------------------------
2104!OLI_CODE
2105!OLI_CODE       perm = 3.0e-08_wp * (minval(phi))**3
2106!OLI_CODE     
2107!OLI_CODE    END SUBROUTINE permeability_phi
2108!OLI_CODE   
2109!OLI_CODE #else
2110!OLI_CODE    !!----------------------------------------------------------------------
2111!OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model
2112!OLI_CODE    !!----------------------------------------------------------------------
2113!OLI_CODE CONTAINS
2114!OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine
2115!OLI_CODE    END SUBROUTINE lim_mp_init   
2116!OLI_CODE    SUBROUTINE lim_mp               ! Empty routine
2117!OLI_CODE    END SUBROUTINE lim_mp
2118!OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine
2119!OLI_CODE    END SUBROUTINE compute_mp_topo       
2120!OLI_CODE    SUBROUTINE pond_area            ! Empty routine
2121!OLI_CODE    END SUBROUTINE pond_area   
2122!OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine
2123!OLI_CODE    END SUBROUTINE calc_hpond   
2124!OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine
2125!OLI_CODE    END SUBROUTINE permeability_phy   
2126!OLI_CODE #endif
2127!OLI_CODE    !!======================================================================
2128!OLI_CODE END MODULE limmp_topo
2129!OLI_CODE
2130#else
2131   !!----------------------------------------------------------------------
2132   !!   Default option          Empty module           NO LIM sea-ice model
2133   !!----------------------------------------------------------------------
2134CONTAINS
2135   SUBROUTINE lim_mp_init     ! Empty routine
2136   END SUBROUTINE lim_mp_init
2137   SUBROUTINE lim_mp          ! Empty routine
2138   END SUBROUTINE lim_mp     
2139   SUBROUTINE lim_mp_topo     ! Empty routine
2140   END SUBROUTINE lim_mp_topo
2141   SUBROUTINE lim_mp_cesm     ! Empty routine
2142   END SUBROUTINE lim_mp_cesm
2143   SUBROUTINE lim_mp_area     ! Empty routine
2144   END SUBROUTINE lim_mp_area
2145   SUBROUTINE lim_mp_perm     ! Empty routine
2146   END SUBROUTINE lim_mp_perm
2147#endif 
2148
2149   !!======================================================================
2150END MODULE limmp 
Note: See TracBrowser for help on using the repository browser.