New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icethd_pnd.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/icethd_pnd.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago

Caught up with trunk rev 14020...

  • Property svn:keywords set to Id
File size: 65.5 KB
Line 
1MODULE icethd_pnd
2   !!======================================================================
3   !!                     ***  MODULE  icethd_pnd   ***
4   !!   sea-ice: Melt ponds on top of sea ice
5   !!======================================================================
6   !! history :       !  2012     (O. Lecomte)       Adaptation from Flocco and Turner
7   !!                 !  2017     (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation
8   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3' :                                     SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_thd_pnd_init : some initialization and namelist read
15   !!   ice_thd_pnd      : main calling routine
16   !!----------------------------------------------------------------------
17   USE phycst         ! physical constants
18   USE dom_oce        ! ocean space and time domain
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE icetab         ! sea-ice: 1D <==> 2D transformation
22   USE sbc_ice        ! surface energy budget
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_thd_pnd_init    ! routine called by icestp.F90
34   PUBLIC   ice_thd_pnd         ! routine called by icestp.F90
35
36   INTEGER ::              nice_pnd    ! choice of the type of pond scheme
37   !                                   ! associated indices:
38   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme
39   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme
40   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme
41   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme
42
43   !--------------------------------------------------------------------------
44   ! Diagnostics for pond volume per area
45   !
46   ! dV/dt = mlt + drn + lid + rnf
47   ! mlt   = input from surface melting
48   ! drn   = drainage through brine network
49   ! lid   = lid growth & melt
50   ! rnf   = runoff (water directly removed out of surface melting + overflow)
51   !
52   ! In topo mode, the pond water lost because it is in the snow is not included in the budget
53   ! In level mode, all terms are incorporated
54   !
55   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_mlt       ! meltwater pond volume input      [kg/m2/s]
56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_drn       ! pond volume lost by drainage     [-]
57   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_lid       ! exchange with lid / refreezing   [-]
58   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_rnf       ! meltwater pond lost to runoff    [-]
59   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_mlt_1d    ! meltwater pond volume input      [-]
60   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_drn_1d    ! pond volume lost by drainage     [-]
61   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_lid_1d    ! exchange with lid / refreezing   [-]
62   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_rnf_1d    ! meltwater pond lost to runoff    [-]
63
64   !! * Substitutions
65#  include "do_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
68   !! $Id$
69   !! Software governed by the CeCILL license (see ./LICENSE)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE ice_thd_pnd
74
75      !!-------------------------------------------------------------------
76      !!               ***  ROUTINE ice_thd_pnd   ***
77      !!
78      !! ** Purpose :   change melt pond fraction and thickness
79      !!
80      !! ** Note    : Melt ponds affect only radiative transfer for now
81      !!              No heat, no salt.
82      !!              The current diagnostics lacks a contribution from drainage
83      !!-------------------------------------------------------------------
84      INTEGER ::   ji, jj, jl        ! loop indices
85      !!-------------------------------------------------------------------
86
87      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) )
88      ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) )
89      !
90      diag_dvpn_mlt (:,:) = 0._wp   ;   diag_dvpn_drn (:,:) = 0._wp
91      diag_dvpn_lid (:,:) = 0._wp   ;   diag_dvpn_rnf (:,:) = 0._wp
92      diag_dvpn_mlt_1d(:) = 0._wp   ;   diag_dvpn_drn_1d(:) = 0._wp
93      diag_dvpn_lid_1d(:) = 0._wp   ;   diag_dvpn_rnf_1d(:) = 0._wp
94
95      !-------------------------------------
96      !  Remove ponds where ice has vanished
97      !-------------------------------------
98      at_i(:,:) = SUM( a_i, dim=3 )
99      !
100      DO jl = 1, jpl
101         DO_2D( 1, 1, 1, 1 )
102            IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN
103               wfx_pnd  (ji,jj)    = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice
104               a_ip     (ji,jj,jl) = 0._wp
105               v_ip     (ji,jj,jl) = 0._wp
106               v_il     (ji,jj,jl) = 0._wp
107               h_ip     (ji,jj,jl) = 0._wp
108               h_il     (ji,jj,jl) = 0._wp
109               a_ip_frac(ji,jj,jl) = 0._wp
110            ENDIF
111         END_2D
112      END DO
113
114      !------------------------------
115      !  Identify grid cells with ice
116      !------------------------------
117      npti = 0   ;   nptidx(:) = 0
118      DO_2D( 1, 1, 1, 1 )
119         IF( at_i(ji,jj) >= epsi10 ) THEN
120            npti = npti + 1
121            nptidx( npti ) = (jj - 1) * jpi + ji
122         ENDIF
123      END_2D
124
125      !------------------------------------
126      !  Select melt pond scheme to be used
127      !------------------------------------
128      IF( npti > 0 ) THEN
129         SELECT CASE ( nice_pnd )
130            !
131         CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==!
132            !
133         CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==!
134            !
135         CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==!
136            !
137         END SELECT
138      ENDIF
139
140      !------------------------------------
141      !  Diagnostics
142      !------------------------------------
143      CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting
144      CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid
145      CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage
146      CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow
147      !
148      DEALLOCATE( diag_dvpn_mlt   , diag_dvpn_lid   , diag_dvpn_drn   , diag_dvpn_rnf    )
149      DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d )
150
151   END SUBROUTINE ice_thd_pnd
152
153
154   SUBROUTINE pnd_CST
155      !!-------------------------------------------------------------------
156      !!                ***  ROUTINE pnd_CST  ***
157      !!
158      !! ** Purpose :   Compute melt pond evolution
159      !!
160      !! ** Method  :   Melt pond fraction and thickness are prescribed
161      !!                to non-zero values when t_su = 0C
162      !!
163      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
164      !!
165      !! ** Note   : Coupling with such melt ponds is only radiative
166      !!             Advection, ridging, rafting... are bypassed
167      !!
168      !! ** References : Bush, G.W., and Trump, D.J. (2017)
169      !!-------------------------------------------------------------------
170      INTEGER  ::   ji, jl    ! loop indices
171      REAL(wp) ::   zdv_pnd   ! Amount of water going into the ponds & lids
172      !!-------------------------------------------------------------------
173      DO jl = 1, jpl
174
175         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d    (1:npti), a_i    (:,:,jl) )
176         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d   (1:npti), t_su   (:,:,jl) )
177         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) )
178         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) )
179         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) )
180         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    )
181
182         DO ji = 1, npti
183            !
184            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji)
185            !
186            IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN
187               h_ip_1d(ji)      = rn_hpnd
188               a_ip_1d(ji)      = rn_apnd * a_i_1d(ji)
189               h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever
190            ELSE
191               h_ip_1d(ji)      = 0._wp
192               a_ip_1d(ji)      = 0._wp
193               h_il_1d(ji)      = 0._wp
194            ENDIF
195            !
196            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
197            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
198            !
199            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd
200            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice
201            !
202         END DO
203
204         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) )
205         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) )
206         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) )
207         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d   (1:npti), v_ip   (:,:,jl) )
208         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d   (1:npti), v_il   (:,:,jl) )
209         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    )
210
211      END DO
212      !
213   END SUBROUTINE pnd_CST
214
215
216   SUBROUTINE pnd_LEV
217      !!-------------------------------------------------------------------
218      !!                ***  ROUTINE pnd_LEV  ***
219      !!
220      !! ** Purpose : Compute melt pond evolution
221      !!
222      !! ** Method  : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
223      !!              We  work with volumes and then redistribute changes into thickness and concentration
224      !!              assuming linear relationship between the two.
225      !!
226      !! ** Action  : - pond growth:      Vp = Vp + dVmelt                                          --- from Holland et al 2012 ---
227      !!                                     dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
228      !!                                        dh_i  = meltwater from ice surface melt
229      !!                                        dh_s  = meltwater from snow melt
230      !!                                        (1-r) = fraction of melt water that is not flushed
231      !!
232      !!              - limtations:       a_ip must not exceed (1-r)*a_i
233      !!                                  h_ip must not exceed 0.5*h_i
234      !!
235      !!              - pond shrinking:
236      !!                       if lids:   Vp = Vp -dH * a_ip
237      !!                                     dH = lid thickness change. Retrieved from this eq.:    --- from Flocco et al 2010 ---
238      !!
239      !!                                                                   rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
240      !!                                                                      H = lid thickness
241      !!                                                                      Lf = latent heat of fusion
242      !!                                                                      Tp = -2C
243      !!
244      !!                                                                And solved implicitely as:
245      !!                                                                   H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
246      !!
247      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 ---
248      !!
249      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi * flush         --- from Flocco et al 2007 ---
250      !!                                     perm = permability of sea-ice                              + correction from Hunke et al 2012 (flush)
251      !!                                     visc = water viscosity
252      !!                                     Hp   = height of top of the pond above sea-level
253      !!                                     Hi   = ice thickness thru which there is flushing
254      !!                                     flush= correction otherwise flushing is excessive
255      !!
256      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness
257      !!
258      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
259      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect
260      !!
261      !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush
262      !!
263      !! ** Note       :   Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds.
264      !!
265      !! ** References :   Flocco and Feltham (JGR, 2007)
266      !!                   Flocco et al       (JGR, 2010)
267      !!                   Holland et al      (J. Clim, 2012)
268      !!                   Hunke et al        (OM 2012)
269      !!-------------------------------------------------------------------
270      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array
271      !!
272      REAL(wp), PARAMETER ::   zaspect =  0.8_wp      ! pond aspect ratio
273      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature
274      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity
275      !!
276      REAL(wp) ::   zfr_mlt, zdv_mlt, zdv_avail       ! fraction and volume of available meltwater retained for melt ponding
277      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes
278      REAL(wp) ::   zdv_pnd                           ! Amount of water going into the ponds & lids
279      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh
280      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed
281      REAL(wp) ::   zdT                               ! zTp-t_su
282      REAL(wp) ::   zsbr, ztmelts                     ! Brine salinity
283      REAL(wp) ::   zperm                             ! permeability of sea ice
284      REAL(wp) ::   zfac, zdum                        ! temporary arrays
285      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
286      !!
287      INTEGER  ::   ji, jk, jl                        ! loop indices
288      !!-------------------------------------------------------------------
289      z1_rhow   = 1._wp / rhow
290      z1_aspect = 1._wp / zaspect
291      z1_Tp     = 1._wp / zTp
292
293      CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d          (1:npti), at_i          )
294      CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d       (1:npti), wfx_pnd       )
295
296      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt )
297      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn )
298      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid )
299      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf )
300
301      DO jl = 1, jpl
302
303         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i (:,:,jl) )
304         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d  (1:npti), h_i (:,:,jl) )
305         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) )
306         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) )
307         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) )
308         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) )
309
310         CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) )
311         CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) )
312
313         DO jk = 1, nlay_i
314            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) )
315            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) )
316         END DO
317
318         !-----------------------
319         ! Melt pond calculations
320         !-----------------------
321         DO ji = 1, npti
322            !
323            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji)
324            !                                                            !----------------------------------------------------!
325            IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN   ! Case ice thickness < rn_himin or tiny ice fraction !
326               !                                                         !----------------------------------------------------!
327               !--- Remove ponds on thin ice or tiny ice fractions
328               a_ip_1d(ji) = 0._wp
329               h_ip_1d(ji) = 0._wp
330               h_il_1d(ji) = 0._wp
331               !                                                         !--------------------------------!
332            ELSE                                                         ! Case ice thickness >= rn_himin !
333               !                                                         !--------------------------------!
334               v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
335               v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
336               !
337               !------------------!
338               ! case ice melting !
339               !------------------!
340               !
341               !--- available meltwater for melt ponding (zdv_avail) ---!
342               zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0
343               zfr_mlt   = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed
344               zdv_mlt   = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors?
345               !
346               !--- overflow ---!
347               !
348               ! area driven overflow
349               !    If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
350               !       a_ip_max = zfr_mlt * a_i
351               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
352               zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
353               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
354
355               ! depth driven overflow
356               !    If pond depth exceeds half the ice thickness then reduce the pond volume
357               !       h_ip_max = 0.5 * h_i
358               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
359               zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
360               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
361
362               !--- Pond growing ---!
363               v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
364               !
365               !--- Lid melting ---!
366               IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
367               !
368               !-------------------!
369               ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
370               !-------------------!
371               !
372               zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
373               !
374               !--- Pond contraction (due to refreezing) ---!
375               IF( ln_pnd_lids ) THEN
376                  !
377                  !--- Lid growing and subsequent pond shrinking ---!
378                  zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
379                     &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
380
381                  ! Lid growing
382                  v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz )
383
384                  ! Pond shrinking
385                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz )
386
387               ELSE
388                  zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp )  ! Holland 2012 (eq. 6)
389                  ! Pond shrinking
390                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz )
391               ENDIF
392               !
393               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
394               ! v_ip     = h_ip * a_ip
395               ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
396               a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
397               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
398               !
399
400               !------------------------------------------------!
401               ! Pond drainage through brine network (flushing) !
402               !------------------------------------------------!
403               ! height of top of the pond above sea-level
404               zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0
405
406               ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
407               DO jk = 1, nlay_i
408                  ! MV Assur is inconsistent with SI3
409                  !!zsbr = - 1.2_wp                                  &
410                  !!   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
411                  !!   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
412                  !!   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3
413                  !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr
414                  ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
415                  ztmelts  = -rTmlt * sz_i_1d(ji,jk)
416                  ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 )
417               END DO
418               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
419
420               ! Do the drainage using Darcy's law
421               zdv_flush   = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012)
422               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0
423               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
424
425               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
426               a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
427               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
428
429               !--- Corrections and lid thickness ---!
430               IF( ln_pnd_lids ) THEN
431                  !--- retrieve lid thickness from volume ---!
432                  IF( a_ip_1d(ji) > 0.01_wp ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
433                  ELSE                               ;   h_il_1d(ji) = 0._wp
434                  ENDIF
435                  !--- remove ponds if lids are much larger than ponds ---!
436                  IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
437                     a_ip_1d(ji) = 0._wp
438                     h_ip_1d(ji) = 0._wp
439                     h_il_1d(ji) = 0._wp
440                  ENDIF
441               ENDIF
442
443               ! diagnostics: dvpnd = mlt+rnf+lid+drn
444               diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow *   zdv_avail             * r1_Dt_ice   ! > 0, surface melt input
445               diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice   ! < 0, runoff
446               diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow *   zdv_frz               * r1_Dt_ice   ! < 0, shrinking
447               diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow *   zdv_flush             * r1_Dt_ice   ! < 0, drainage
448               !
449            ENDIF
450            !
451            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
452            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
453            !
454            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd
455            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice
456            !
457         END DO
458
459         !--------------------------------------------------------------------
460         ! Retrieve 2D arrays
461         !--------------------------------------------------------------------
462         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) )
463         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) )
464         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) )
465         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) )
466         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) )
467         !
468      END DO
469      !
470      CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd )
471      !
472      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        )
473      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        )
474      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        )
475      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )
476      !
477   END SUBROUTINE pnd_LEV
478
479
480
481   SUBROUTINE pnd_TOPO
482
483      !!-------------------------------------------------------------------
484      !!                ***  ROUTINE pnd_TOPO  ***
485      !!
486      !! ** Purpose :   Compute melt pond evolution based on the ice
487      !!                topography inferred from the ice thickness distribution
488      !!
489      !! ** Method  :   This code is initially based on Flocco and Feltham
490      !!                (2007) and Flocco et al. (2010).
491      !!
492      !!                - Calculate available pond water base on surface meltwater
493      !!                - Redistribute water as a function of topography, drain water
494      !!                - Exchange water with the lid
495      !!
496      !! ** Tunable parameters :
497      !!
498      !! ** Note :
499      !!
500      !! ** References
501      !!
502      !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond
503      !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:
504      !!    10.1029/2006JC003836.
505      !!
506      !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of
507      !!    a physically based melt pond scheme into the sea ice component of a
508      !!    climate model.  J. Geophys. Res. 115, C08012,
509      !!    doi: 10.1029/2009JC005568.
510      !!
511      !!-------------------------------------------------------------------
512      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds
513         zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C)
514         zvp_min = 1.e-4_wp          ! minimum pond volume (m)
515
516
517      ! local variables
518      REAL(wp) :: &
519         zdHui,   &      ! change in thickness of ice lid (m)
520         zomega,  &      ! conduction
521         zdTice,  &      ! temperature difference across ice lid (C)
522         zdvice,  &      ! change in ice volume (m)
523         zTavg,   &      ! mean surface temperature across categories (C)
524         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2)
525         zTp,     &      ! pond freezing temperature (C)
526         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3)
527         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding
528         z1_rhow, &      ! inverse water density
529         zv_pnd  , &     ! volume of meltwater contributing to ponds
530         zv_mlt          ! total amount of meltwater produced
531
532      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, &     !! total melt pond water available before redistribution and drainage
533                                        zvolp_res    !! remaining melt pond water available after drainage
534
535      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i
536
537      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
538
539      INTEGER  ::   i_test
540
541      ! Note
542      ! equivalent for CICE translation
543      ! a_ip      -> apond
544      ! a_ip_frac -> apnd
545
546      CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' )
547
548      !---------------------------------------------------------------
549      ! Initialise
550      !---------------------------------------------------------------
551
552      ! Parameters & constants (move to parameters)
553      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3)
554      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline)
555      z1_rhow   = 1._wp / rhow
556
557      ! Set required ice variables (hard-coded here for now)
558!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?)
559
560      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction
561      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area
562      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area
563      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area
564
565      ! thickness
566      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
567      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
568      END WHERE
569      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
570
571      !---------------------------------------------------------------
572      ! Change 2D to 1D
573      !---------------------------------------------------------------
574      ! MV
575      ! a less computing-intensive version would have 2D-1D passage here
576      ! use what we have in iceitd.F90 (incremental remapping)
577
578      !--------------------------------------------------------------
579      ! Collect total available pond water volume
580      !--------------------------------------------------------------
581      ! Assuming that meltwater (+rain in principle) runsoff the surface
582      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction
583      ! I cite her words, they are very talkative
584      ! "grid cells with very little ice cover (and hence more open water area)
585      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water."
586      ! "This results in the same runoff fraction r for each ice category within a grid cell"
587
588      zvolp(:,:) = 0._wp
589
590      DO jl = 1, jpl
591         DO_2D( 1, 1, 1, 1 )
592
593               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
594
595                  !--- Available and contributing meltwater for melt ponding ---!
596                  zv_mlt  = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) &        ! available volume of surface melt water per grid area
597                     &    * z1_rhow * a_i(ji,jj,jl)
598                      ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area
599                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds
600                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl
601
602                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice                     ! diags
603                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice
604
605                  !--- Create possible new ponds
606                  ! if pond does not exist, create new pond over full ice area
607                  !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
608                  IF ( a_ip(ji,jj,jl) < epsi10 ) THEN
609                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl)
610                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE)
611                  ENDIF
612
613                  !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage
614                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness
615                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
616
617                  !--- Total available pond water volume (pre-existing + newly produced)j
618                  zvolp(ji,jj)   = zvolp(ji,jj)   + v_ip(ji,jj,jl)
619!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now
620
621               ENDIF ! a_i
622
623         END_2D
624      END DO ! ji
625
626      !--------------------------------------------------------------
627      ! Redistribute and drain water from ponds
628      !--------------------------------------------------------------
629      CALL ice_thd_pnd_area( zvolp, zvolp_res )
630
631      !--------------------------------------------------------------
632      ! Melt pond lid growth and melt
633      !--------------------------------------------------------------
634
635      IF( ln_pnd_lids ) THEN
636
637         DO_2D( 1, 1, 1, 1 )
638
639            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
640
641               !--------------------------
642               ! Pond lid growth and melt
643               !--------------------------
644               ! Mean surface temperature
645               zTavg = 0._wp
646               DO jl = 1, jpl
647                  zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
648               END DO
649               zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
650
651               DO jl = 1, jpl-1
652
653                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN
654
655                     !----------------------------------------------------------------
656                     ! Lid melting: floating upper ice layer melts in whole or part
657                     !----------------------------------------------------------------
658                     ! Use Tsfc for each category
659                     IF ( t_su(ji,jj,jl) > zTp ) THEN
660
661                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
662
663                        IF ( zdvice > epsi10 ) THEN
664
665                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
666                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -
667                                                                        ! as it is already counted in surface melt
668!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar)
669!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw)
670
671                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
672                           ! ice lid melted and category is pond covered
673                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl)
674!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl)
675                              v_il(ji,jj,jl)   = 0._wp
676                           ENDIF
677                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
678
679                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag
680
681                        ENDIF
682
683                     !----------------------------------------------------------------
684                     ! Freeze pre-existing lid
685                     !----------------------------------------------------------------
686
687                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
688
689                        ! differential growth of base of surface floating ice layer
690                        zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0
691                        zomega = rcnd_i * zdTice / zrhoi_L
692                        zdHui  = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
693                               - v_il(ji,jj,jl) / a_i(ji,jj,jl)
694                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
695
696                        IF ( zdvice > epsi10 ) THEN
697                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
698                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
699!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice
700!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice
701                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
702
703                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag
704
705                        ENDIF
706
707                     ENDIF ! Tsfcn(i,j,n)
708
709                  !----------------------------------------------------------------
710                  ! Freeze new lids
711                  !----------------------------------------------------------------
712                  !  upper ice layer begins to form
713                  ! note: albedo does not change
714
715                  ELSE ! v_il < epsi10
716
717                     ! thickness of newly formed ice
718                     ! the surface temperature of a meltpond is the same as that
719                     ! of the ice underneath (0C), and the thermodynamic surface
720                     ! flux is the same
721
722                     !!! we need net surface energy flux, excluding conduction
723                     !!! fsurf is summed over categories in CICE
724                     !!! we have the category-dependent flux, let us use it ?
725                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)
726                     zdHui  = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp )
727                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
728                     IF ( zdvice > epsi10 ) THEN
729                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
730                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
731
732                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag
733!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice
734!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice
735                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar
736                     ENDIF
737
738                  ENDIF  ! v_il
739
740               END DO ! jl
741
742            ELSE  ! remove ponds on thin ice
743
744               v_ip(ji,jj,:) = 0._wp
745               v_il(ji,jj,:) = 0._wp
746!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
747!                 zvolp(ji,jj)    = 0._wp
748
749            ENDIF
750
751         END_2D
752
753      ENDIF ! ln_pnd_lids
754
755      !---------------------------------------------------------------
756      ! Clean-up variables (probably duplicates what icevar would do)
757      !---------------------------------------------------------------
758      ! MV comment
759      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il
760      ! icevar should recompute all other variables (if needed at all)
761
762      DO jl = 1, jpl
763
764         DO_2D( 1, 1, 1, 1 )
765
766!              ! zap lids on small ponds
767!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
768!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN
769!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small
770!              ENDIF
771
772               ! recalculate equivalent pond variables
773               IF ( a_ip(ji,jj,jl) > epsi10) THEN
774                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_i(ji,jj,jl)
775                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
776                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
777               ENDIF
778!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar
779!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar
780!              ENDIF
781
782         END_2D
783
784      END DO   ! jl
785
786
787   END SUBROUTINE pnd_TOPO
788
789
790   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp )
791
792       !!-------------------------------------------------------------------
793       !!                ***  ROUTINE ice_thd_pnd_area ***
794       !!
795       !! ** Purpose : Given the total volume of available pond water,
796       !!              redistribute and drain water
797       !!
798       !! ** Method
799       !!
800       !-----------|
801       !           |
802       !           |-----------|
803       !___________|___________|______________________________________sea-level
804       !           |           |
805       !           |           |---^--------|
806       !           |           |   |        |
807       !           |           |   |        |-----------|              |-------
808       !           |           |   | alfan  |           |              |
809       !           |           |   |        |           |--------------|
810       !           |           |   |        |           |              |
811       !---------------------------v-------------------------------------------
812       !           |           |   ^        |           |              |
813       !           |           |   |        |           |--------------|
814       !           |           |   | betan  |           |              |
815       !           |           |   |        |-----------|              |-------
816       !           |           |   |        |
817       !           |           |---v------- |
818       !           |           |
819       !           |-----------|
820       !           |
821       !-----------|
822       !
823       !!
824       !!------------------------------------------------------------------
825
826       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
827          zvolp,                                       &  ! total available pond water
828          zdvolp                                          ! remaining meltwater after redistribution
829
830       INTEGER ::  &
831          ns,      &
832          m_index, &
833          permflag
834
835       REAL (wp), DIMENSION(jpl) :: &
836          hicen, &
837          hsnon, &
838          asnon, &
839          alfan, &
840          betan, &
841          cum_max_vol, &
842          reduced_aicen
843
844       REAL (wp), DIMENSION(0:jpl) :: &
845          cum_max_vol_tmp
846
847       REAL (wp) :: &
848          hpond, &
849          drain, &
850          floe_weight, &
851          pressure_head, &
852          hsl_rel, &
853          deltah, &
854          perm, &
855          msno
856
857       REAL (wp), parameter :: &
858          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s
859
860      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds
861         zvp_min = 1.e-4_wp          ! minimum pond volume (m)
862
863      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
864
865       a_ip(:,:,:) = 0._wp
866       h_ip(:,:,:) = 0._wp
867
868       DO_2D( 1, 1, 1, 1 )
869
870             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
871
872        !-------------------------------------------------------------------
873        ! initialize
874        !-------------------------------------------------------------------
875
876        DO jl = 1, jpl
877
878           !----------------------------------------
879           ! compute the effective snow fraction
880           !----------------------------------------
881
882           IF (a_i(ji,jj,jl) < epsi10)  THEN
883              hicen(jl) =  0._wp
884              hsnon(jl) =  0._wp
885              reduced_aicen(jl) = 0._wp
886              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
887           ELSE
888              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
889              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl)
890              reduced_aicen(jl) = 1._wp ! n=jpl
891
892              !js: initial code in NEMO_DEV
893              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) &
894              !                     * (-0.024_wp*hicen(jl) + 0.832_wp)
895
896              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
897              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &
898                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp))
899
900              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical)
901              ! MV should check whether this makes sense to have the same effective snow fraction in here
902              ! OLI: it probably doesn't
903           END IF
904
905  ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
906  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
907  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
908  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
909  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
910  ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
911  ! alfan = 60% of the ice volume) in each category lies above the reference line,
912  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
913
914  ! MV:
915  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE
916  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
917
918  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
919
920           alfan(jl) = 0.6 * hicen(jl)
921           betan(jl) = 0.4 * hicen(jl)
922
923           cum_max_vol(jl)     = 0._wp
924           cum_max_vol_tmp(jl) = 0._wp
925
926        END DO ! jpl
927
928        cum_max_vol_tmp(0) = 0._wp
929        drain = 0._wp
930        zdvolp(ji,jj) = 0._wp
931
932        !----------------------------------------------------------
933        ! Drain overflow water, update pond fraction and volume
934        !----------------------------------------------------------
935
936        !--------------------------------------------------------------------------
937        ! the maximum amount of water that can be contained up to each ice category
938        !--------------------------------------------------------------------------
939        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
940        ! Then the excess volume cum_max_vol(jl) drains out of the system
941        ! It should be added to wfx_pnd_out
942
943        DO jl = 1, jpl-1 ! last category can not hold any volume
944
945           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN
946
947              ! total volume in level including snow
948              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + &
949                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl))
950
951              ! subtract snow solid volumes from lower categories in current level
952              DO ns = 1, jl
953                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
954                    - rhos/rhow * &     ! free air fraction that can be filled by water
955                      asnon(ns)  * &    ! effective areal fraction of snow in that category
956                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
957              END DO
958
959           ELSE ! assume higher categories unoccupied
960              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1)
961           END IF
962           !IF (cum_max_vol_tmp(jl) < z0) THEN
963           !   CALL abort_ice('negative melt pond volume')
964           !END IF
965        END DO
966        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
967        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
968
969        !----------------------------------------------------------------
970        ! is there more meltwater than can be held in the floe?
971        !----------------------------------------------------------------
972        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN
973           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
974           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
975
976           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice)
977
978           zdvolp(ji,jj) = drain         ! this is the drained water
979           IF (zvolp(ji,jj) < epsi10) THEN
980              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
981              zvolp(ji,jj) = 0._wp
982           END IF
983        END IF
984
985        ! height and area corresponding to the remaining volume
986        ! routine leaves zvolp unchanged
987        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
988
989        DO jl = 1, m_index
990           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update
991           !                                         !  volume instead, no ?
992           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2
993           a_ip(ji,jj,jl) = reduced_aicen(jl)
994           ! in practise, pond fraction depends on the empirical snow fraction
995           ! so in turn on ice thickness
996        END DO
997        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
998
999        !------------------------------------------------------------------------
1000        ! Drainage through brine network (permeability)
1001        !------------------------------------------------------------------------
1002        !!! drainage due to ice permeability - Darcy's law
1003
1004        ! sea water level
1005        msno = 0._wp
1006        DO jl = 1 , jpl
1007          msno = msno + v_s(ji,jj,jl) * rhos
1008        END DO
1009        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj)
1010        hsl_rel = floe_weight / rho0 &
1011                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) )
1012
1013        deltah = hpond - hsl_rel
1014        pressure_head = grav * rho0 * max(deltah, 0._wp)
1015
1016        ! drain if ice is permeable
1017        permflag = 0
1018
1019        IF (pressure_head > 0._wp) THEN
1020           DO jl = 1, jpl-1
1021              IF ( hicen(jl) /= 0._wp ) THEN
1022
1023              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2
1024
1025                 perm = 0._wp ! MV ugly dummy patch
1026                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof
1027                 IF (perm > 0._wp) permflag = 1
1028
1029                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / &
1030                                          (viscosity*hicen(jl))
1031                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
1032                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
1033
1034                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded)
1035
1036                 IF (zvolp(ji,jj) < epsi10) THEN
1037                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
1038                    zvolp(ji,jj) = 0._wp
1039                 END IF
1040             END IF
1041          END DO
1042
1043           ! adjust melt pond dimensions
1044           IF (permflag > 0) THEN
1045              ! recompute pond depth
1046             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
1047              DO jl = 1, m_index
1048                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1)
1049                 a_ip(ji,jj,jl) = reduced_aicen(jl)
1050              END DO
1051              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
1052           END IF
1053        END IF ! pressure_head
1054
1055        !-------------------------------
1056        ! remove water from the snow
1057        !-------------------------------
1058        !------------------------------------------------------------------------
1059        ! total melt pond volume in category does not include snow volume
1060        ! snow in melt ponds is not melted
1061        !------------------------------------------------------------------------
1062
1063        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell
1064        ! how much, so I did not diagnose it
1065        ! so if there is a problem here, nobody is going to see it...
1066
1067
1068        ! Calculate pond volume for lower categories
1069        DO jl = 1,m_index-1
1070           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
1071                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
1072        END DO
1073
1074        ! Calculate pond volume for highest category = remaining pond volume
1075
1076        ! The following is completely unclear to Martin at least
1077        ! Could we redefine properly and recode in a more readable way ?
1078
1079        ! m_index = last category with melt pond
1080
1081        IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water
1082
1083        IF (m_index > 1) THEN
1084          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN
1085             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1))
1086          ELSE
1087             v_ip(ji,jj,m_index) = 0._wp
1088             h_ip(ji,jj,m_index) = 0._wp
1089             a_ip(ji,jj,m_index) = 0._wp
1090             ! If remaining pond volume is negative reduce pond volume of
1091             ! lower category
1092             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) &
1093              v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj)
1094          END IF
1095        END IF
1096
1097        DO jl = 1,m_index
1098           IF (a_ip(ji,jj,jl) > epsi10) THEN
1099               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
1100           ELSE
1101              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl)
1102              h_ip(ji,jj,jl) = 0._wp
1103              v_ip(ji,jj,jl)  = 0._wp
1104              a_ip(ji,jj,jl) = 0._wp
1105           END IF
1106        END DO
1107        DO jl = m_index+1, jpl
1108           h_ip(ji,jj,jl) = 0._wp
1109           a_ip(ji,jj,jl) = 0._wp
1110           v_ip(ji,jj,jl) = 0._wp
1111        END DO
1112
1113           ENDIF
1114
1115     END_2D
1116
1117    END SUBROUTINE ice_thd_pnd_area
1118
1119
1120    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
1121       !!-------------------------------------------------------------------
1122       !!                ***  ROUTINE ice_thd_pnd_depth  ***
1123       !!
1124       !! ** Purpose :   Compute melt pond depth
1125       !!-------------------------------------------------------------------
1126
1127       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1128          aicen, &
1129          asnon, &
1130          hsnon, &
1131          alfan, &
1132          cum_max_vol
1133
1134       REAL (wp), INTENT(IN) :: &
1135          zvolp
1136
1137       REAL (wp), INTENT(OUT) :: &
1138          hpond
1139
1140       INTEGER, INTENT(OUT) :: &
1141          m_index
1142
1143       INTEGER :: n, ns
1144
1145       REAL (wp), DIMENSION(0:jpl+1) :: &
1146          hitl, &
1147          aicetl
1148
1149       REAL (wp) :: &
1150          rem_vol, &
1151          area, &
1152          vol, &
1153          tmp, &
1154          z0   = 0.0_wp
1155
1156       !----------------------------------------------------------------
1157       ! hpond is zero if zvolp is zero - have we fully drained?
1158       !----------------------------------------------------------------
1159
1160       IF (zvolp < epsi10) THEN
1161        hpond = z0
1162        m_index = 0
1163       ELSE
1164
1165        !----------------------------------------------------------------
1166        ! Calculate the category where water fills up to
1167        !----------------------------------------------------------------
1168
1169        !----------|
1170        !          |
1171        !          |
1172        !          |----------|                                     -- --
1173        !__________|__________|_________________________________________ ^
1174        !          |          |             rem_vol     ^                | Semi-filled
1175        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
1176        !          |          |          |              |
1177        !          |          |          |              |hpond
1178        !          |          |          |----------|   |     |-------
1179        !          |          |          |          |   |     |
1180        !          |          |          |          |---v-----|
1181        !          |          | m_index  |          |         |
1182        !-------------------------------------------------------------
1183
1184        m_index = 0  ! 1:m_index categories have water in them
1185        DO n = 1, jpl
1186           IF (zvolp <= cum_max_vol(n)) THEN
1187              m_index = n
1188              IF (n == 1) THEN
1189                 rem_vol = zvolp
1190              ELSE
1191                 rem_vol = zvolp - cum_max_vol(n-1)
1192              END IF
1193              exit ! to break out of the loop
1194           END IF
1195        END DO
1196        m_index = min(jpl-1, m_index)
1197
1198        !----------------------------------------------------------------
1199        ! semi-filled layer may have m_index different snow in it
1200        !----------------------------------------------------------------
1201
1202        !-----------------------------------------------------------  ^
1203        !                                                             |  alfan(m_index+1)
1204        !                                                             |
1205        !hitl(3)-->                             |----------|          |
1206        !hitl(2)-->                |------------| * * * * *|          |
1207        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
1208        !hitl(0)-->-------------------------------------------------  |  ^
1209        !                various snow from lower categories          |  |alfa(m_index)
1210
1211        ! hitl - heights of the snow layers from thinner and current categories
1212        ! aicetl - area of each snow depth in this layer
1213
1214        hitl(:) = z0
1215        aicetl(:) = z0
1216        DO n = 1, m_index
1217           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1218                                  alfan(m_index+1) - alfan(m_index)), z0)
1219           aicetl(n) = asnon(n)
1220
1221           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1222        END DO
1223
1224        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1225        aicetl(m_index+1) = z0
1226
1227        !----------------------------------------------------------------
1228        ! reorder array according to hitl
1229        ! snow heights not necessarily in height order
1230        !----------------------------------------------------------------
1231
1232        DO ns = 1, m_index+1
1233           DO n = 0, m_index - ns + 1
1234              IF (hitl(n) > hitl(n+1)) THEN ! swap order
1235                 tmp = hitl(n)
1236                 hitl(n) = hitl(n+1)
1237                 hitl(n+1) = tmp
1238                 tmp = aicetl(n)
1239                 aicetl(n) = aicetl(n+1)
1240                 aicetl(n+1) = tmp
1241              END IF
1242           END DO
1243        END DO
1244
1245        !----------------------------------------------------------------
1246        ! divide semi-filled layer into set of sublayers each vertically homogenous
1247        !----------------------------------------------------------------
1248
1249        !hitl(3)----------------------------------------------------------------
1250        !                                                       | * * * * * * * *
1251        !                                                       |* * * * * * * * *
1252        !hitl(2)----------------------------------------------------------------
1253        !                                    | * * * * * * * *  | * * * * * * * *
1254        !                                    |* * * * * * * * * |* * * * * * * * *
1255        !hitl(1)----------------------------------------------------------------
1256        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *
1257        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
1258        !hitl(0)----------------------------------------------------------------
1259        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)
1260
1261        ! move up over layers incrementing volume
1262        DO n = 1, m_index+1
1263
1264           area = sum(aicetl(:)) - &                 ! total area of sub-layer
1265                (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1266
1267           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
1268
1269           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
1270              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
1271
1272              exit
1273           ELSE  ! still in sub-layer below the sub-layer with the depth
1274              rem_vol = rem_vol - vol
1275           END IF
1276
1277        END DO
1278
1279       END IF
1280
1281    END SUBROUTINE ice_thd_pnd_depth
1282
1283
1284    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
1285       !!-------------------------------------------------------------------
1286       !!                ***  ROUTINE ice_thd_pnd_perm ***
1287       !!
1288       !! ** Purpose :   Determine the liquid fraction of brine in the ice
1289       !!                and its permeability
1290       !!-------------------------------------------------------------------
1291
1292       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1293          ticen,  &     ! internal ice temperature (K)
1294          salin         ! salinity (ppt)     !js: ppt according to cice
1295
1296       REAL (wp), INTENT(OUT) :: &
1297          perm      ! permeability
1298
1299       REAL (wp) ::   &
1300          Sbr       ! brine salinity
1301
1302       REAL (wp), DIMENSION(nlay_i) ::   &
1303          Tin, &    ! ice temperature
1304          phi       ! liquid fraction
1305
1306       INTEGER :: k
1307
1308       !-----------------------------------------------------------------
1309       ! Compute ice temperatures from enthalpies using quadratic formula
1310       !-----------------------------------------------------------------
1311
1312       DO k = 1,nlay_i
1313          Tin(k) = ticen(k) - rt0   !js: from K to degC
1314       END DO
1315
1316       !-----------------------------------------------------------------
1317       ! brine salinity and liquid fraction
1318       !-----------------------------------------------------------------
1319
1320       DO k = 1, nlay_i
1321
1322          Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
1323          ! Best expression to date is that one (Vancoppenolle et al JGR 2019)
1324          ! Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3
1325          phi(k) = salin(k) / Sbr
1326
1327       END DO
1328
1329       !-----------------------------------------------------------------
1330       ! permeability
1331       !-----------------------------------------------------------------
1332
1333       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
1334
1335   END SUBROUTINE ice_thd_pnd_perm
1336
1337   SUBROUTINE ice_thd_pnd_init
1338      !!-------------------------------------------------------------------
1339      !!                  ***  ROUTINE ice_thd_pnd_init   ***
1340      !!
1341      !! ** Purpose : Physical constants and parameters linked to melt ponds
1342      !!              over sea ice
1343      !!
1344      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond
1345      !!               parameter values called at the first timestep (nit000)
1346      !!
1347      !! ** input   :   Namelist namthd_pnd
1348      !!-------------------------------------------------------------------
1349      INTEGER  ::   ios, ioptio   ! Local integer
1350      !!
1351      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, &
1352         &                          ln_pnd_CST , rn_apnd, rn_hpnd,         &
1353         &                          ln_pnd_TOPO,                           &
1354         &                          ln_pnd_lids, ln_pnd_alb
1355      !!-------------------------------------------------------------------
1356      !
1357      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
1358901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
1359      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
1360902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
1361      IF(lwm) WRITE ( numoni, namthd_pnd )
1362      !
1363      IF(lwp) THEN                        ! control print
1364         WRITE(numout,*)
1365         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
1366         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
1367         WRITE(numout,*) '   Namelist namicethd_pnd:'
1368         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
1369         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO
1370         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV
1371         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
1372         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
1373         WRITE(numout,*) '            Pond flushing efficiency                              rn_pnd_flush = ', rn_pnd_flush
1374         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST
1375         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
1376         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
1377         WRITE(numout,*) '         Frozen lids on top of melt ponds                         ln_pnd_lids  = ', ln_pnd_lids
1378         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
1379      ENDIF
1380      !
1381      !                             !== set the choice of ice pond scheme ==!
1382      ioptio = 0
1383      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
1384      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
1385      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF
1386      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF
1387      IF( ioptio /= 1 )   &
1388         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' )
1389      !
1390      SELECT CASE( nice_pnd )
1391      CASE( np_pndNO )
1392         IF( ln_pnd_alb  ) THEN ; ln_pnd_alb  = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' )  ; ENDIF
1393         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
1394      CASE( np_pndCST )
1395         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
1396      END SELECT
1397      !
1398   END SUBROUTINE ice_thd_pnd_init
1399
1400#else
1401   !!----------------------------------------------------------------------
1402   !!   Default option          Empty module          NO SI3 sea-ice model
1403   !!----------------------------------------------------------------------
1404#endif
1405
1406   !!======================================================================
1407END MODULE icethd_pnd
Note: See TracBrowser for help on using the repository browser.