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

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

impact of melt ponds on albedo, first run

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