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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icethd_pnd.F90

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

trunk: solve ticket #2625

  • Property svn:keywords set to Id
File size: 65.8 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( nn_hls, nn_hls, nn_hls, nn_hls )
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( nn_hls, nn_hls, nn_hls, nn_hls )
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     = 273._wp       , & ! temperature difference for freeze-up (K)
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_ini , &   !! total melt pond water available before redistribution and drainage
533                                        zvolp     , &   !! total melt pond water volume
534                                        zvolp_res       !! remaining melt pond water available after drainage
535
536      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i
537
538      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
539
540      INTEGER  ::   i_test
541
542      ! Note
543      ! equivalent for CICE translation
544      ! a_ip      -> apond
545      ! a_ip_frac -> apnd
546
547      CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' )
548
549      !---------------------------------------------------------------
550      ! Initialise
551      !---------------------------------------------------------------
552
553      ! Parameters & constants (move to parameters)
554      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3)
555      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline)
556      z1_rhow   = 1._wp / rhow
557
558      ! Set required ice variables (hard-coded here for now)
559!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?)
560
561      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction
562      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area
563      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area
564      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area
565
566      ! thickness
567      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
568      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
569      END WHERE
570      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
571
572      !---------------------------------------------------------------
573      ! Change 2D to 1D
574      !---------------------------------------------------------------
575      ! MV
576      ! a less computing-intensive version would have 2D-1D passage here
577      ! use what we have in iceitd.F90 (incremental remapping)
578
579      !--------------------------------------------------------------
580      ! Collect total available pond water volume
581      !--------------------------------------------------------------
582      ! Assuming that meltwater (+rain in principle) runsoff the surface
583      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction
584      ! I cite her words, they are very talkative
585      ! "grid cells with very little ice cover (and hence more open water area)
586      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water."
587      ! "This results in the same runoff fraction r for each ice category within a grid cell"
588
589      zvolp(:,:) = 0._wp
590
591      DO jl = 1, jpl
592         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
593
594               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
595
596                  !--- Available and contributing meltwater for melt ponding ---!
597                  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
598                     &    * z1_rhow * a_i(ji,jj,jl)
599                      ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area
600                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds
601                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl
602
603                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice                     ! diags
604                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice
605
606                  !--- Create possible new ponds
607                  ! if pond does not exist, create new pond over full ice area
608                  !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
609                  IF ( a_ip(ji,jj,jl) < epsi10 ) THEN
610                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl)
611                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE)
612                  ENDIF
613
614                  !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage
615                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness
616                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
617
618                  !--- Total available pond water volume (pre-existing + newly produced)j
619                  zvolp(ji,jj)   = zvolp(ji,jj)   + v_ip(ji,jj,jl)
620!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now
621
622               ENDIF ! a_i
623
624         END_2D
625      END DO ! ji
626
627      zvolp_ini(:,:) = zvolp(:,:)
628
629      !--------------------------------------------------------------
630      ! Redistribute and drain water from ponds
631      !--------------------------------------------------------------
632      CALL ice_thd_pnd_area( zvolp, zvolp_res )
633
634      !--------------------------------------------------------------
635      ! Melt pond lid growth and melt
636      !--------------------------------------------------------------
637
638      IF( ln_pnd_lids ) THEN
639
640         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
641
642            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
643
644               !--------------------------
645               ! Pond lid growth and melt
646               !--------------------------
647               ! Mean surface temperature
648               zTavg = 0._wp
649               DO jl = 1, jpl
650                  zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
651               END DO
652               zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
653
654               DO jl = 1, jpl-1
655
656                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN
657
658                     !----------------------------------------------------------------
659                     ! Lid melting: floating upper ice layer melts in whole or part
660                     !----------------------------------------------------------------
661                     ! Use Tsfc for each category
662                     IF ( t_su(ji,jj,jl) > zTp ) THEN
663
664                        zdvice = MIN( -dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
665
666                        IF ( zdvice > epsi10 ) THEN
667
668                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
669                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -
670                                                                        ! as it is already counted in surface melt
671!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar)
672!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw)
673
674                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
675                           ! ice lid melted and category is pond covered
676                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl)
677!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl)
678                              v_il(ji,jj,jl)   = 0._wp
679                           ENDIF
680                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
681
682                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag
683
684                        ENDIF
685
686                     !----------------------------------------------------------------
687                     ! Freeze pre-existing lid
688                     !----------------------------------------------------------------
689
690                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
691
692                        ! differential growth of base of surface floating ice layer
693                        zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0
694                        zomega = rcnd_i * zdTice / zrhoi_L
695                        zdHui  = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
696                               - v_il(ji,jj,jl) / a_i(ji,jj,jl)
697                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
698
699                        IF ( zdvice > epsi10 ) THEN
700                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
701                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
702!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice
703!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice
704                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
705
706                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag
707
708                        ENDIF
709
710                     ENDIF ! Tsfcn(i,j,n)
711
712                  !----------------------------------------------------------------
713                  ! Freeze new lids
714                  !----------------------------------------------------------------
715                  !  upper ice layer begins to form
716                  ! note: albedo does not change
717
718                  ELSE ! v_il < epsi10
719
720                     ! thickness of newly formed ice
721                     ! the surface temperature of a meltpond is the same as that
722                     ! of the ice underneath (0C), and the thermodynamic surface
723                     ! flux is the same
724
725                     !!! we need net surface energy flux, excluding conduction
726                     !!! fsurf is summed over categories in CICE
727                     !!! we have the category-dependent flux, let us use it ?
728                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)
729                     zdHui  = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp )
730                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
731                     IF ( zdvice > epsi10 ) THEN
732                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
733                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
734
735                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag
736!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice
737!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice
738                        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
739                     ENDIF
740
741                  ENDIF  ! v_il
742
743               END DO ! jl
744
745            ELSE  ! remove ponds on thin ice
746
747               v_ip(ji,jj,:) = 0._wp
748               v_il(ji,jj,:) = 0._wp
749!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
750!                 zvolp(ji,jj)    = 0._wp
751
752            ENDIF
753
754         END_2D
755
756      ENDIF ! ln_pnd_lids
757
758      !---------------------------------------------------------------
759      ! Clean-up variables (probably duplicates what icevar would do)
760      !---------------------------------------------------------------
761      ! MV comment
762      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il
763      ! icevar should recompute all other variables (if needed at all)
764
765      DO jl = 1, jpl
766
767         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
768
769!              ! zap lids on small ponds
770!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
771!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN
772!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small
773!              ENDIF
774
775               ! recalculate equivalent pond variables
776               IF ( a_ip(ji,jj,jl) > epsi10) THEN
777                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
778                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
779                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
780               ENDIF
781!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar
782!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar
783!              ENDIF
784
785         END_2D
786
787      END DO   ! jl
788
789
790   END SUBROUTINE pnd_TOPO
791
792
793   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp )
794
795       !!-------------------------------------------------------------------
796       !!                ***  ROUTINE ice_thd_pnd_area ***
797       !!
798       !! ** Purpose : Given the total volume of available pond water,
799       !!              redistribute and drain water
800       !!
801       !! ** Method
802       !!
803       !-----------|
804       !           |
805       !           |-----------|
806       !___________|___________|______________________________________sea-level
807       !           |           |
808       !           |           |---^--------|
809       !           |           |   |        |
810       !           |           |   |        |-----------|              |-------
811       !           |           |   | alfan  |           |              |
812       !           |           |   |        |           |--------------|
813       !           |           |   |        |           |              |
814       !---------------------------v-------------------------------------------
815       !           |           |   ^        |           |              |
816       !           |           |   |        |           |--------------|
817       !           |           |   | betan  |           |              |
818       !           |           |   |        |-----------|              |-------
819       !           |           |   |        |
820       !           |           |---v------- |
821       !           |           |
822       !           |-----------|
823       !           |
824       !-----------|
825       !
826       !!
827       !!------------------------------------------------------------------
828
829       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
830          zvolp,                                       &  ! total available pond water
831          zdvolp                                          ! remaining meltwater after redistribution
832
833       INTEGER ::  &
834          ns,      &
835          m_index, &
836          permflag
837
838       REAL (wp), DIMENSION(jpl) :: &
839          hicen, &
840          hsnon, &
841          asnon, &
842          alfan, &
843          betan, &
844          cum_max_vol, &
845          reduced_aicen
846
847       REAL (wp), DIMENSION(0:jpl) :: &
848          cum_max_vol_tmp
849
850       REAL (wp) :: &
851          hpond, &
852          drain, &
853          floe_weight, &
854          pressure_head, &
855          hsl_rel, &
856          deltah, &
857          perm, &
858          msno
859
860       REAL (wp), parameter :: &
861          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s
862
863      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds
864         zvp_min = 1.e-4_wp          ! minimum pond volume (m)
865
866      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
867
868       a_ip(:,:,:) = 0._wp
869       h_ip(:,:,:) = 0._wp
870
871       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
872
873             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
874
875        !-------------------------------------------------------------------
876        ! initialize
877        !-------------------------------------------------------------------
878
879        DO jl = 1, jpl
880
881           !----------------------------------------
882           ! compute the effective snow fraction
883           !----------------------------------------
884
885           IF (a_i(ji,jj,jl) < epsi10)  THEN
886              hicen(jl) =  0._wp
887              hsnon(jl) =  0._wp
888              reduced_aicen(jl) = 0._wp
889              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
890           ELSE
891              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
892              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl)
893              reduced_aicen(jl) = 1._wp ! n=jpl
894
895              !js: initial code in NEMO_DEV
896              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) &
897              !                     * (-0.024_wp*hicen(jl) + 0.832_wp)
898
899              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
900              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &
901                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp))
902
903              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical)
904              ! MV should check whether this makes sense to have the same effective snow fraction in here
905              ! OLI: it probably doesn't
906           END IF
907
908  ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
909  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
910  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
911  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
912  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
913  ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
914  ! alfan = 60% of the ice volume) in each category lies above the reference line,
915  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
916
917  ! MV:
918  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE
919  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
920
921  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
922
923           alfan(jl) = 0.6 * hicen(jl)
924           betan(jl) = 0.4 * hicen(jl)
925
926           cum_max_vol(jl)     = 0._wp
927           cum_max_vol_tmp(jl) = 0._wp
928
929        END DO ! jpl
930
931        cum_max_vol_tmp(0) = 0._wp
932        drain = 0._wp
933        zdvolp(ji,jj) = 0._wp
934
935        !----------------------------------------------------------
936        ! Drain overflow water, update pond fraction and volume
937        !----------------------------------------------------------
938
939        !--------------------------------------------------------------------------
940        ! the maximum amount of water that can be contained up to each ice category
941        !--------------------------------------------------------------------------
942        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
943        ! Then the excess volume cum_max_vol(jl) drains out of the system
944        ! It should be added to wfx_pnd_out
945
946        DO jl = 1, jpl-1 ! last category can not hold any volume
947
948           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN
949
950              ! total volume in level including snow
951              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + &
952                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl))
953
954              ! subtract snow solid volumes from lower categories in current level
955              DO ns = 1, jl
956                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
957                    - rhos/rhow * &     ! free air fraction that can be filled by water
958                      asnon(ns)  * &    ! effective areal fraction of snow in that category
959                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
960              END DO
961
962           ELSE ! assume higher categories unoccupied
963              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1)
964           END IF
965           !IF (cum_max_vol_tmp(jl) < z0) THEN
966           !   CALL abort_ice('negative melt pond volume')
967           !END IF
968        END DO
969        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
970        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
971
972        !----------------------------------------------------------------
973        ! is there more meltwater than can be held in the floe?
974        !----------------------------------------------------------------
975        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN
976           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
977           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
978
979           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice)
980
981           zdvolp(ji,jj) = drain         ! this is the drained water
982           IF (zvolp(ji,jj) < epsi10) THEN
983              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
984              zvolp(ji,jj) = 0._wp
985           END IF
986        END IF
987
988        ! height and area corresponding to the remaining volume
989        ! routine leaves zvolp unchanged
990        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
991
992        DO jl = 1, m_index
993           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update
994           !                                         !  volume instead, no ?
995           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2
996           a_ip(ji,jj,jl) = reduced_aicen(jl)
997           ! in practise, pond fraction depends on the empirical snow fraction
998           ! so in turn on ice thickness
999        END DO
1000        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
1001
1002        !------------------------------------------------------------------------
1003        ! Drainage through brine network (permeability)
1004        !------------------------------------------------------------------------
1005        !!! drainage due to ice permeability - Darcy's law
1006
1007        ! sea water level
1008        msno = 0._wp
1009        DO jl = 1 , jpl
1010          msno = msno + v_s(ji,jj,jl) * rhos
1011        END DO
1012        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj)
1013        hsl_rel = floe_weight / rho0 &
1014                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) )
1015
1016        deltah = hpond - hsl_rel
1017        pressure_head = grav * rho0 * max(deltah, 0._wp)
1018
1019        ! drain if ice is permeable
1020        permflag = 0
1021
1022        IF (pressure_head > 0._wp) THEN
1023           DO jl = 1, jpl-1
1024              IF ( hicen(jl) /= 0._wp ) THEN
1025
1026              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2
1027
1028                 perm = 0._wp ! MV ugly dummy patch
1029                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof
1030                 IF (perm > 0._wp) permflag = 1
1031
1032                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / &
1033                                          (viscosity*hicen(jl))
1034                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
1035                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
1036
1037                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded)
1038
1039                 IF (zvolp(ji,jj) < epsi10) THEN
1040                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
1041                    zvolp(ji,jj) = 0._wp
1042                 END IF
1043             END IF
1044          END DO
1045
1046           ! adjust melt pond dimensions
1047           IF (permflag > 0) THEN
1048              ! recompute pond depth
1049             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
1050              DO jl = 1, m_index
1051                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1)
1052                 a_ip(ji,jj,jl) = reduced_aicen(jl)
1053              END DO
1054              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
1055           END IF
1056        END IF ! pressure_head
1057
1058        !-------------------------------
1059        ! remove water from the snow
1060        !-------------------------------
1061        !------------------------------------------------------------------------
1062        ! total melt pond volume in category does not include snow volume
1063        ! snow in melt ponds is not melted
1064        !------------------------------------------------------------------------
1065
1066        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell
1067        ! how much, so I did not diagnose it
1068        ! so if there is a problem here, nobody is going to see it...
1069
1070
1071        ! Calculate pond volume for lower categories
1072        DO jl = 1,m_index-1
1073           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
1074                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
1075        END DO
1076
1077        ! Calculate pond volume for highest category = remaining pond volume
1078
1079        ! The following is completely unclear to Martin at least
1080        ! Could we redefine properly and recode in a more readable way ?
1081
1082        ! m_index = last category with melt pond
1083
1084        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
1085
1086        IF (m_index > 1) THEN
1087          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN
1088             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1))
1089          ELSE
1090             v_ip(ji,jj,m_index) = 0._wp
1091             h_ip(ji,jj,m_index) = 0._wp
1092             a_ip(ji,jj,m_index) = 0._wp
1093             ! If remaining pond volume is negative reduce pond volume of
1094             ! lower category
1095             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) &
1096              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)
1097          END IF
1098        END IF
1099
1100        DO jl = 1,m_index
1101           IF (a_ip(ji,jj,jl) > epsi10) THEN
1102               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
1103           ELSE
1104              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl)
1105              h_ip(ji,jj,jl) = 0._wp
1106              v_ip(ji,jj,jl)  = 0._wp
1107              a_ip(ji,jj,jl) = 0._wp
1108           END IF
1109        END DO
1110        DO jl = m_index+1, jpl
1111           h_ip(ji,jj,jl) = 0._wp
1112           a_ip(ji,jj,jl) = 0._wp
1113           v_ip(ji,jj,jl) = 0._wp
1114        END DO
1115
1116           ENDIF
1117
1118     END_2D
1119
1120    END SUBROUTINE ice_thd_pnd_area
1121
1122
1123    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
1124       !!-------------------------------------------------------------------
1125       !!                ***  ROUTINE ice_thd_pnd_depth  ***
1126       !!
1127       !! ** Purpose :   Compute melt pond depth
1128       !!-------------------------------------------------------------------
1129
1130       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
1131          aicen, &
1132          asnon, &
1133          hsnon, &
1134          alfan, &
1135          cum_max_vol
1136
1137       REAL (wp), INTENT(IN) :: &
1138          zvolp
1139
1140       REAL (wp), INTENT(OUT) :: &
1141          hpond
1142
1143       INTEGER, INTENT(OUT) :: &
1144          m_index
1145
1146       INTEGER :: n, ns
1147
1148       REAL (wp), DIMENSION(0:jpl+1) :: &
1149          hitl, &
1150          aicetl
1151
1152       REAL (wp) :: &
1153          rem_vol, &
1154          area, &
1155          vol, &
1156          tmp, &
1157          z0   = 0.0_wp
1158
1159       !----------------------------------------------------------------
1160       ! hpond is zero if zvolp is zero - have we fully drained?
1161       !----------------------------------------------------------------
1162
1163       IF (zvolp < epsi10) THEN
1164        hpond = z0
1165        m_index = 0
1166       ELSE
1167
1168        !----------------------------------------------------------------
1169        ! Calculate the category where water fills up to
1170        !----------------------------------------------------------------
1171
1172        !----------|
1173        !          |
1174        !          |
1175        !          |----------|                                     -- --
1176        !__________|__________|_________________________________________ ^
1177        !          |          |             rem_vol     ^                | Semi-filled
1178        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
1179        !          |          |          |              |
1180        !          |          |          |              |hpond
1181        !          |          |          |----------|   |     |-------
1182        !          |          |          |          |   |     |
1183        !          |          |          |          |---v-----|
1184        !          |          | m_index  |          |         |
1185        !-------------------------------------------------------------
1186
1187        m_index = 0  ! 1:m_index categories have water in them
1188        DO n = 1, jpl
1189           IF (zvolp <= cum_max_vol(n)) THEN
1190              m_index = n
1191              IF (n == 1) THEN
1192                 rem_vol = zvolp
1193              ELSE
1194                 rem_vol = zvolp - cum_max_vol(n-1)
1195              END IF
1196              exit ! to break out of the loop
1197           END IF
1198        END DO
1199        m_index = min(jpl-1, m_index)
1200
1201        !----------------------------------------------------------------
1202        ! semi-filled layer may have m_index different snow in it
1203        !----------------------------------------------------------------
1204
1205        !-----------------------------------------------------------  ^
1206        !                                                             |  alfan(m_index+1)
1207        !                                                             |
1208        !hitl(3)-->                             |----------|          |
1209        !hitl(2)-->                |------------| * * * * *|          |
1210        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
1211        !hitl(0)-->-------------------------------------------------  |  ^
1212        !                various snow from lower categories          |  |alfa(m_index)
1213
1214        ! hitl - heights of the snow layers from thinner and current categories
1215        ! aicetl - area of each snow depth in this layer
1216
1217        hitl(:) = z0
1218        aicetl(:) = z0
1219        DO n = 1, m_index
1220           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
1221                                  alfan(m_index+1) - alfan(m_index)), z0)
1222           aicetl(n) = asnon(n)
1223
1224           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
1225        END DO
1226
1227        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
1228        aicetl(m_index+1) = z0
1229
1230        !----------------------------------------------------------------
1231        ! reorder array according to hitl
1232        ! snow heights not necessarily in height order
1233        !----------------------------------------------------------------
1234
1235        DO ns = 1, m_index+1
1236           DO n = 0, m_index - ns + 1
1237              IF (hitl(n) > hitl(n+1)) THEN ! swap order
1238                 tmp = hitl(n)
1239                 hitl(n) = hitl(n+1)
1240                 hitl(n+1) = tmp
1241                 tmp = aicetl(n)
1242                 aicetl(n) = aicetl(n+1)
1243                 aicetl(n+1) = tmp
1244              END IF
1245           END DO
1246        END DO
1247
1248        !----------------------------------------------------------------
1249        ! divide semi-filled layer into set of sublayers each vertically homogenous
1250        !----------------------------------------------------------------
1251
1252        !hitl(3)----------------------------------------------------------------
1253        !                                                       | * * * * * * * *
1254        !                                                       |* * * * * * * * *
1255        !hitl(2)----------------------------------------------------------------
1256        !                                    | * * * * * * * *  | * * * * * * * *
1257        !                                    |* * * * * * * * * |* * * * * * * * *
1258        !hitl(1)----------------------------------------------------------------
1259        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *
1260        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
1261        !hitl(0)----------------------------------------------------------------
1262        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)
1263
1264        ! move up over layers incrementing volume
1265        DO n = 1, m_index+1
1266
1267           area = sum(aicetl(:)) - &                 ! total area of sub-layer
1268                (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
1269
1270           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
1271
1272           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
1273              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
1274
1275              exit
1276           ELSE  ! still in sub-layer below the sub-layer with the depth
1277              rem_vol = rem_vol - vol
1278           END IF
1279
1280        END DO
1281
1282       END IF
1283
1284    END SUBROUTINE ice_thd_pnd_depth
1285
1286
1287    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
1288       !!-------------------------------------------------------------------
1289       !!                ***  ROUTINE ice_thd_pnd_perm ***
1290       !!
1291       !! ** Purpose :   Determine the liquid fraction of brine in the ice
1292       !!                and its permeability
1293       !!-------------------------------------------------------------------
1294
1295       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
1296          ticen,  &     ! internal ice temperature (K)
1297          salin         ! salinity (ppt)     !js: ppt according to cice
1298
1299       REAL (wp), INTENT(OUT) :: &
1300          perm      ! permeability
1301
1302       REAL (wp) ::   &
1303          Sbr       ! brine salinity
1304
1305       REAL (wp), DIMENSION(nlay_i) ::   &
1306          Tin, &    ! ice temperature
1307          phi       ! liquid fraction
1308
1309       INTEGER :: k
1310
1311       !-----------------------------------------------------------------
1312       ! Compute ice temperatures from enthalpies using quadratic formula
1313       !-----------------------------------------------------------------
1314
1315       DO k = 1,nlay_i
1316          Tin(k) = ticen(k) - rt0   !js: from K to degC
1317       END DO
1318
1319       !-----------------------------------------------------------------
1320       ! brine salinity and liquid fraction
1321       !-----------------------------------------------------------------
1322
1323       DO k = 1, nlay_i
1324
1325          Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
1326          ! Best expression to date is that one (Vancoppenolle et al JGR 2019)
1327          ! Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3
1328          phi(k) = salin(k) / Sbr
1329
1330       END DO
1331
1332       !-----------------------------------------------------------------
1333       ! permeability
1334       !-----------------------------------------------------------------
1335
1336       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
1337
1338   END SUBROUTINE ice_thd_pnd_perm
1339
1340   SUBROUTINE ice_thd_pnd_init
1341      !!-------------------------------------------------------------------
1342      !!                  ***  ROUTINE ice_thd_pnd_init   ***
1343      !!
1344      !! ** Purpose : Physical constants and parameters linked to melt ponds
1345      !!              over sea ice
1346      !!
1347      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond
1348      !!               parameter values called at the first timestep (nit000)
1349      !!
1350      !! ** input   :   Namelist namthd_pnd
1351      !!-------------------------------------------------------------------
1352      INTEGER  ::   ios, ioptio   ! Local integer
1353      !!
1354      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, &
1355         &                          ln_pnd_CST , rn_apnd, rn_hpnd,         &
1356         &                          ln_pnd_TOPO,                           &
1357         &                          ln_pnd_lids, ln_pnd_alb
1358      !!-------------------------------------------------------------------
1359      !
1360      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
1361901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
1362      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
1363902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
1364      IF(lwm) WRITE ( numoni, namthd_pnd )
1365      !
1366      IF(lwp) THEN                        ! control print
1367         WRITE(numout,*)
1368         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
1369         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
1370         WRITE(numout,*) '   Namelist namicethd_pnd:'
1371         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
1372         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO
1373         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV
1374         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
1375         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
1376         WRITE(numout,*) '            Pond flushing efficiency                              rn_pnd_flush = ', rn_pnd_flush
1377         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST
1378         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
1379         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
1380         WRITE(numout,*) '         Frozen lids on top of melt ponds                         ln_pnd_lids  = ', ln_pnd_lids
1381         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
1382      ENDIF
1383      !
1384      !                             !== set the choice of ice pond scheme ==!
1385      ioptio = 0
1386      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
1387      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
1388      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF
1389      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF
1390      IF( ioptio /= 1 )   &
1391         & 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)' )
1392      !
1393      SELECT CASE( nice_pnd )
1394      CASE( np_pndNO )
1395         IF( ln_pnd_alb  ) THEN ; ln_pnd_alb  = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' )  ; ENDIF
1396         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
1397      CASE( np_pndCST )
1398         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
1399      END SELECT
1400      !
1401   END SUBROUTINE ice_thd_pnd_init
1402
1403#else
1404   !!----------------------------------------------------------------------
1405   !!   Default option          Empty module          NO SI3 sea-ice model
1406   !!----------------------------------------------------------------------
1407#endif
1408
1409   !!======================================================================
1410END MODULE icethd_pnd
Note: See TracBrowser for help on using the repository browser.