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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_pnd.F90

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

4.0-HEAD: change flushing param in SI3 to match better with melt pond observations

  • Property svn:keywords set to Id
File size: 21.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   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_pnd_init    ! routine called by icestp.F90
32   PUBLIC   ice_thd_pnd         ! routine called by icestp.F90
33
34   INTEGER ::              nice_pnd    ! choice of the type of pond scheme
35   !                                   ! associated indices:
36   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
37   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant ice pond scheme
38   INTEGER, PARAMETER ::   np_pndLEV = 2   ! Level ice pond scheme
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE ice_thd_pnd
50      !!-------------------------------------------------------------------
51      !!               ***  ROUTINE ice_thd_pnd   ***
52      !!               
53      !! ** Purpose :   change melt pond fraction and thickness
54      !!               
55      !!-------------------------------------------------------------------
56      !
57      SELECT CASE ( nice_pnd )
58      !
59      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
60         !
61      CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==!
62         !
63      END SELECT
64      !
65   END SUBROUTINE ice_thd_pnd 
66
67
68   SUBROUTINE pnd_CST 
69      !!-------------------------------------------------------------------
70      !!                ***  ROUTINE pnd_CST  ***
71      !!
72      !! ** Purpose :   Compute melt pond evolution
73      !!
74      !! ** Method  :   Melt pond fraction and thickness are prescribed
75      !!                to non-zero values when t_su = 0C
76      !!
77      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
78      !!               
79      !! ** Note   : Coupling with such melt ponds is only radiative
80      !!             Advection, ridging, rafting... are bypassed
81      !!
82      !! ** References : Bush, G.W., and Trump, D.J. (2017)
83      !!-------------------------------------------------------------------
84      INTEGER  ::   ji        ! loop indices
85      REAL(wp) ::   zdv_pnd   ! Amount of water going into the ponds & lids
86      !!-------------------------------------------------------------------
87      DO ji = 1, npti
88         !
89         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji)
90         !
91         IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN
92            h_ip_1d(ji)      = rn_hpnd   
93            a_ip_1d(ji)      = rn_apnd * a_i_1d(ji)
94            h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever
95         ELSE
96            h_ip_1d(ji)      = 0._wp   
97            a_ip_1d(ji)      = 0._wp
98            h_il_1d(ji)      = 0._wp
99         ENDIF
100         !
101         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd
102         wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice
103         !
104      END DO
105      !
106   END SUBROUTINE pnd_CST
107
108
109   SUBROUTINE pnd_LEV
110      !!-------------------------------------------------------------------
111      !!                ***  ROUTINE pnd_LEV  ***
112      !!
113      !! ** Purpose : Compute melt pond evolution
114      !!
115      !! ** Method  : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
116      !!              We  work with volumes and then redistribute changes into thickness and concentration
117      !!              assuming linear relationship between the two.
118      !!
119      !! ** Action  : - pond growth:      Vp = Vp + dVmelt                                          --- from Holland et al 2012 ---
120      !!                                     dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
121      !!                                        dh_i  = meltwater from ice surface melt
122      !!                                        dh_s  = meltwater from snow melt
123      !!                                        (1-r) = fraction of melt water that is not flushed
124      !!
125      !!              - limtations:       a_ip must not exceed (1-r)*a_i
126      !!                                  h_ip must not exceed 0.5*h_i
127      !!
128      !!              - pond shrinking:
129      !!                       if lids:   Vp = Vp -dH * a_ip
130      !!                                     dH = lid thickness change. Retrieved from this eq.:    --- from Flocco et al 2010 ---
131      !!
132      !!                                                                   rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
133      !!                                                                      H = lid thickness
134      !!                                                                      Lf = latent heat of fusion
135      !!                                                                      Tp = -2C
136      !!
137      !!                                                                And solved implicitely as:
138      !!                                                                   H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
139      !!
140      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 ---
141      !!
142      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi * flush         --- from Flocco et al 2007 ---
143      !!                                     perm = permability of sea-ice                              + correction from Hunke et al 2012 (flush)
144      !!                                     visc = water viscosity
145      !!                                     Hp   = height of top of the pond above sea-level
146      !!                                     Hi   = ice thickness thru which there is flushing
147      !!                                     flush= correction otherwise flushing is excessive
148      !!
149      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness
150      !!
151      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
152      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect
153      !!
154      !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
155      !!
156      !! ** Note       :   mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds.
157      !!
158      !! ** References :   Flocco and Feltham (JGR, 2007)
159      !!                   Flocco et al       (JGR, 2010)
160      !!                   Holland et al      (J. Clim, 2012)
161      !!                   Hunke et al        (OM 2012)
162      !!-------------------------------------------------------------------
163      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array
164      !!
165      REAL(wp), PARAMETER ::   zaspect =  0.8_wp      ! pond aspect ratio
166      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature
167      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity
168      REAL(wp), PARAMETER ::   zflush  =  0.1_wp      ! tuning param to reduce flushing (from Hunke et al 2012)
169      !!
170      REAL(wp) ::   zfr_mlt, zdv_mlt, zdv_avail       ! fraction and volume of available meltwater retained for melt ponding
171      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes
172      REAL(wp) ::   zdv_pnd                           ! Amount of water going into the ponds & lids
173      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh
174      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed
175      REAL(wp) ::   zdT                               ! zTp-t_su
176      REAL(wp) ::   zsbr, ztmelts                     ! Brine salinity
177      REAL(wp) ::   zperm                             ! permeability of sea ice
178      REAL(wp) ::   zfac, zdum                        ! temporary arrays
179      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
180      !!
181      INTEGER  ::   ji, jk                            ! loop indices
182      !!-------------------------------------------------------------------
183      z1_rhow   = 1._wp / rhow 
184      z1_aspect = 1._wp / zaspect
185      z1_Tp     = 1._wp / zTp 
186
187      DO ji = 1, npti
188         !
189         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji)
190         !                                                            !----------------------------------------------------!
191         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 !
192            !                                                         !----------------------------------------------------!
193            !--- Remove ponds on thin ice or tiny ice fractions
194            a_ip_1d(ji) = 0._wp
195            h_ip_1d(ji) = 0._wp
196            h_il_1d(ji) = 0._wp
197            !                                                         !--------------------------------!
198         ELSE                                                         ! Case ice thickness >= rn_himin !
199            !                                                         !--------------------------------!
200            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
201            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
202            !
203            !------------------!
204            ! case ice melting !
205            !------------------!
206            !
207            !--- available meltwater for melt ponding (zdv_avail) ---!
208            zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0
209            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
210            zdv_mlt   = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors?
211            !
212            !--- overflow ---!
213            !
214            ! area driven overflow
215            !    If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
216            !       a_ip_max = zfr_mlt * a_i
217            !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
218            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
219            zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
220
221            ! depth driven overflow
222            !    If pond depth exceeds half the ice thickness then reduce the pond volume
223            !       h_ip_max = 0.5 * h_i
224            !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
225            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
226            zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
227           
228            !--- Pond growing ---!
229            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
230            !
231            !--- Lid melting ---!
232            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
233            !
234            !-------------------!
235            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
236            !-------------------!
237            !
238            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
239            !
240            !--- Pond contraction (due to refreezing) ---!
241            IF( ln_pnd_lids ) THEN
242               !
243               !--- Lid growing and subsequent pond shrinking ---!
244               zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
245                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
246               
247               ! Lid growing
248               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz )
249               
250               ! Pond shrinking
251               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz )
252
253            ELSE
254               zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp )  ! Holland 2012 (eq. 6)
255               ! Pond shrinking
256               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz )
257            ENDIF
258            !
259            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
260            ! v_ip     = h_ip * a_ip
261            ! 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)
262            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
263            h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
264           
265            !------------------------------------------------!           
266            ! Pond drainage through brine network (flushing) !
267            !------------------------------------------------!
268            ! height of top of the pond above sea-level
269            zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
270           
271            ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
272            DO jk = 1, nlay_i
273               ! MV Assur is inconsistent with SI3
274               !!zsbr = - 1.2_wp                                  &
275               !!   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
276               !!   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
277               !!   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3
278               !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr
279               ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
280               ztmelts  = -rTmlt * sz_i_1d(ji,jk)
281               ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 )
282            END DO
283            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
284           
285            ! Do the drainage using Darcy's law
286            zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * zflush ! zflush comes from Hunke et al. (2012)
287            zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0
288            v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
289           
290            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
291            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
292            h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
293
294            !--- Corrections and lid thickness ---!
295            IF( ln_pnd_lids ) THEN
296               !--- retrieve lid thickness from volume ---!
297               IF( a_ip_1d(ji) > 0.01_wp ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
298               ELSE                               ;   h_il_1d(ji) = 0._wp
299               ENDIF
300               !--- remove ponds if lids are much larger than ponds ---!
301               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
302                  a_ip_1d(ji) = 0._wp
303                  h_ip_1d(ji) = 0._wp
304                  h_il_1d(ji) = 0._wp
305               ENDIF
306            ENDIF
307
308!!$            ! diagnostics: dvpnd = mlt+rnf+frz+drn
309!!$            diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow *   zdv_avail             * r1_rdtice   ! > 0, surface melt input
310!!$            diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_rdtice   ! < 0, runoff
311!!$            diag_dvpn_frz_1d(ji) = diag_dvpn_frz_1d(ji) + rhow *   zdv_frz               * r1_rdtice   ! < 0, shrinking
312!!$            diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow *   zdv_flush             * r1_rdtice   ! < 0, drainage
313            !
314         ENDIF
315         !
316         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd
317         wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice
318         !
319      END DO
320      !
321   END SUBROUTINE pnd_LEV
322
323
324   SUBROUTINE ice_thd_pnd_init 
325      !!-------------------------------------------------------------------
326      !!                  ***  ROUTINE ice_thd_pnd_init   ***
327      !!
328      !! ** Purpose : Physical constants and parameters linked to melt ponds
329      !!              over sea ice
330      !!
331      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
332      !!               parameter values called at the first timestep (nit000)
333      !!
334      !! ** input   :   Namelist namthd_pnd 
335      !!-------------------------------------------------------------------
336      INTEGER  ::   ios, ioptio   ! Local integer
337      !!
338      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &
339         &                          ln_pnd_CST , rn_apnd, rn_hpnd,         &
340         &                          ln_pnd_lids, ln_pnd_alb
341      !!-------------------------------------------------------------------
342      !
343      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
344      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
345901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
346      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
347      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
348902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
349      IF(lwm) WRITE ( numoni, namthd_pnd )
350      !
351      IF(lwp) THEN                        ! control print
352         WRITE(numout,*)
353         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
354         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
355         WRITE(numout,*) '   Namelist namicethd_pnd:'
356         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
357         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV
358         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
359         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
360         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST
361         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
362         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
363         WRITE(numout,*) '         Frozen lids on top of melt ponds                         ln_pnd_lids  = ', ln_pnd_lids
364         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
365      ENDIF
366      !
367      !                             !== set the choice of ice pond scheme ==!
368      ioptio = 0
369      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
370      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
371      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF
372      IF( ioptio /= 1 )   &
373         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )
374      !
375      SELECT CASE( nice_pnd )
376      CASE( np_pndNO )         
377         IF( ln_pnd_alb  ) THEN ; ln_pnd_alb  = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' )  ; ENDIF
378         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
379      CASE( np_pndCST )         
380         IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
381      END SELECT
382      !
383   END SUBROUTINE ice_thd_pnd_init
384   
385#else
386   !!----------------------------------------------------------------------
387   !!   Default option          Empty module          NO SI3 sea-ice model
388   !!----------------------------------------------------------------------
389#endif 
390
391   !!======================================================================
392END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.