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

Last change on this file since 8411 was 8411, checked in by clem, 3 years ago

continue changing names

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