source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90 @ 8106

Last change on this file since 8106 was 8106, checked in by vancop, 4 years ago

Follow up on melt ponds and albedo

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