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

Last change on this file since 8125 was 8125, checked in by vancop, 3 years ago

Melt pond podifications

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