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

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

Further melt pond work

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