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

source: NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90 @ 13912

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

Topographic ponds

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