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

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

Melt pond interfaces practically operational

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