source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90 @ 12720

Last change on this file since 12720 was 12720, checked in by clem, 8 months ago

implementation of ice pond lids (before debugging)

  • Property svn:keywords set to Id
File size: 20.9 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 pond scheme
38   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
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
54      !!               
55      !! ** Method  :   brut force
56      !!-------------------------------------------------------------------
57      !
58      SELECT CASE ( nice_pnd )
59      !
60      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
61         !
62      CASE (np_pndH12)   ;   CALL pnd_H12    !==  Holland et al 2012 melt ponds  ==!
63         !
64      END SELECT
65      !
66   END SUBROUTINE ice_thd_pnd 
67
68
69   SUBROUTINE pnd_CST 
70      !!-------------------------------------------------------------------
71      !!                ***  ROUTINE pnd_CST  ***
72      !!
73      !! ** Purpose :   Compute melt pond evolution
74      !!
75      !! ** Method  :   Melt pond fraction and thickness are prescribed
76      !!                to non-zero values when t_su = 0C
77      !!
78      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
79      !!               
80      !! ** Note   : Coupling with such melt ponds is only radiative
81      !!             Advection, ridging, rafting... are bypassed
82      !!
83      !! ** References : Bush, G.W., and Trump, D.J. (2017)
84      !!-------------------------------------------------------------------
85      INTEGER  ::   ji        ! loop indices
86      !!-------------------------------------------------------------------
87      DO ji = 1, npti
88         !
89         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
90            a_ip_frac_1d(ji) = rn_apnd
91            a_ip_eff_1d(ji)  = rn_apnd
92            h_ip_1d(ji)      = rn_hpnd   
93            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
94            h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever
95         ELSE
96            a_ip_frac_1d(ji) = 0._wp
97            a_ip_eff_1d(ji)  = 0._wp
98            h_ip_1d(ji)      = 0._wp   
99            a_ip_1d(ji)      = 0._wp
100            h_il_1d(ji)      = 0._wp
101         ENDIF
102         !
103      END DO
104      !
105   END SUBROUTINE pnd_CST
106
107
108   SUBROUTINE pnd_H12
109      !!-------------------------------------------------------------------
110      !!                ***  ROUTINE pnd_H12  ***
111      !!
112      !! ** Purpose : Compute melt pond evolution
113      !!
114      !! ** Method  : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
115      !!              We  work with volumes and then redistribute changes into thickness and concentration
116      !!              assuming linear relationship between the two.
117      !!
118      !! ** Action  : - pond growth:      Vp = Vp + dVmelt                                          --- from Holland et al 2012 ---
119      !!                                     dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
120      !!                                        dh_i  = meltwater from ice surface melt
121      !!                                        dh_s  = meltwater from snow melt
122      !!                                        (1-r) = fraction of melt water that is not flushed
123      !!
124      !!              - limtations:       a_ip must not exceed (1-r)*a_i
125      !!                                  h_ip must not exceed 0.5*h_i
126      !!
127      !!              - pond shrinking:
128      !!                       if lids:   Vp = Vp -dH * a_ip
129      !!                                     dH = lid thickness change. Retrieved from this eq.:    --- from Flocco et al 2010 ---
130      !!
131      !!                                                                   rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
132      !!                                                                      H = lid thickness
133      !!                                                                      Lf = latent heat of fusion
134      !!                                                                      Tp = -2C
135      !!
136      !!                                                                And solved implicitely as:
137      !!                                                                   H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
138      !!
139      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 ---
140      !!
141      !!              - Overflow:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 ---
142      !!                                     perm = permability of sea-ice
143      !!                                     visc = water viscosity
144      !!                                     Hp   = height of top of the pond above sea-level
145      !!                                     Hi   = ice thickness thru which there is flushing
146      !!
147      !!
148      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness
149      !!
150      !!              - effective pond area: to be used for albedo
151      !!
152      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
153      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect
154      !!
155      !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
156      !!
157      !! ** Note       :   mostly stolen from CICE
158      !!
159      !! ** References :   Flocco and Feltham (JGR, 2007)
160      !!                   Flocco et al       (JGR, 2010)
161      !!                   Holland et al      (J. Clim, 2012)
162      !!-------------------------------------------------------------------
163      REAL(wp), DIMENSION(nlay_i) ::   zperm          ! Permeability of sea ice
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 ::   zhl_max =  0.015_wp    ! pond lid thickness above which the ponds disappear from the albedo calculation
169      REAL(wp), PARAMETER ::   zhl_min =  0.005_wp    ! pond lid thickness below which the full pond area is used in the albedo calculation
170      !!
171      REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding
172      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes
173      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh
174      REAL(wp) ::   v_ip_max                          ! max pond volume allowed
175      REAL(wp) ::   zdT                               ! zTp-t_su
176      REAL(wp) ::   zsbr                              ! Brine salinity
177      REAL(wp) ::   zfac, zdum                        ! temporary arrays
178      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
179      !!
180      INTEGER  ::   ji, jk                            ! loop indices
181      !!-------------------------------------------------------------------
182      z1_rhow   = 1._wp / rhow 
183      z1_aspect = 1._wp / zaspect
184      z1_Tp     = 1._wp / zTp 
185
186      DO ji = 1, npti
187         !                                                            !----------------------------------------------------!
188         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
189            !                                                         !----------------------------------------------------!
190            !--- Remove ponds on thin ice or tiny ice fractions
191            a_ip_1d(ji)      = 0._wp
192            a_ip_frac_1d(ji) = 0._wp
193            h_ip_1d(ji)      = 0._wp
194            h_il_1d(ji)      = 0._wp
195            !
196            ! clem: problem with conservation or not ?
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 ---!
208            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
209            zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) in H12 = fraction of melt water that is not flushed
210            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?
211            !
212            !--- overflow ---!
213            ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
214            !    a_ip_max = zfr_mlt * a_i
215            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
216            v_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
217            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) )
218
219            ! If pond depth exceeds half the ice thickness then reduce the pond volume
220            !    h_ip_max = 0.5 * h_i
221            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
222            v_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
223            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) )
224           
225            !--- Pond growing ---!
226            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
227            !
228            !--- Lid melting ---!
229            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
230            !
231            !--- mass flux ---!
232            IF( zdv_mlt > 0._wp ) THEN
233               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
234               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
235               !
236               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
237               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
238               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
239            ENDIF
240
241            !-------------------!
242            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
243            !-------------------!
244            !
245            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
246            !
247            !--- Pond contraction (due to refreezing) ---!
248            IF( ln_pnd_lids ) THEN
249               !
250               !--- Lid growing and subsequent pond shrinking ---!
251               zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
252                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
253               
254               ! Lid growing
255               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
256               
257               ! Pond shrinking
258               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
259
260            ELSE
261               ! Pond shrinking
262               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
263            ENDIF
264            !
265            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
266            ! v_ip     = h_ip * a_ip
267            ! 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)
268            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) )
269            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
270            h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji)
271
272            !---------------!           
273            ! Pond flushing !
274            !---------------!
275            IF( ln_pnd_flush ) THEN
276               ! height of top of the pond above sea-level
277               zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_frac_1d(ji) ) ) * r1_rau0
278
279               ! Calculate the permeability of the ice (Assur 1958)
280               DO jk = 1, nlay_i
281                  zsbr = - 1.2_wp                                  &
282                     &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
283                     &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
284                     &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010)
285                  zperm(jk) = MAX( 0._wp, 3.e-08_wp * (sz_i_1d(ji,jk) / zsbr)**3 )
286               END DO
287
288               ! Do the drainage using Darcy's law
289               zdv_flush   = -MINVAL(zperm(:)) * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
290               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
291               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
292               
293               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
294               a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) )
295               a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
296               h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji)
297
298            ENDIF
299
300            !--- Corrections and lid thickness ---!
301            IF( ln_pnd_lids ) THEN
302               !--- remove ponds if lids are much larger than ponds ---!
303               IF ( v_il_1d(ji) > v_ip_1d(ji) * 10._wp ) THEN
304                  a_ip_1d(ji)      = 0._wp
305                  a_ip_frac_1d(ji) = 0._wp
306                  h_ip_1d(ji)      = 0._wp
307                  v_il_1d(ji)      = 0._wp
308               ENDIF
309               !--- retrieve lid thickness from volume ---!
310               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
311               ELSE                              ;   h_il_1d(ji) = 0._wp ;
312               ENDIF
313            ENDIF
314            !
315         ENDIF
316         
317      END DO
318
319      !-------------------------------------------------!           
320      ! How much melt pond is exposed to the atmosphere !
321      !-------------------------------------------------!           
322      ! Calculate the melt pond effective area (used for albedo)
323      WHERE    ( h_il_1d(1:npti) <= zhl_min )   ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti)       ! lid is very thin.  Expose all the pond
324      ELSEWHERE( h_il_1d(1:npti) >= zhl_max )   ;   a_ip_eff_1d(1:npti) = 0._wp                      ! lid is very thick. Cover all the pond up with ice and snow
325      ELSEWHERE                                 ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti) * &   ! lid is in between. Expose part of the pond
326         &                                                                ( h_il_1d(1:npti) - zhl_min ) / ( zhl_max - zhl_min )
327      END WHERE
328      !
329   END SUBROUTINE pnd_H12
330
331
332   SUBROUTINE ice_thd_pnd_init 
333      !!-------------------------------------------------------------------
334      !!                  ***  ROUTINE ice_thd_pnd_init   ***
335      !!
336      !! ** Purpose : Physical constants and parameters linked to melt ponds
337      !!              over sea ice
338      !!
339      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
340      !!               parameter values called at the first timestep (nit000)
341      !!
342      !! ** input   :   Namelist namthd_pnd 
343      !!-------------------------------------------------------------------
344      INTEGER  ::   ios, ioptio   ! Local integer
345      !!
346      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, &
347         &                          ln_pnd_CST, rn_apnd, rn_hpnd, &
348         &                          ln_pnd_alb
349      !!-------------------------------------------------------------------
350      !
351      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
352      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
353901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
354      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
355      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
356902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
357      IF(lwm) WRITE ( numoni, namthd_pnd )
358      !
359      IF(lwp) THEN                        ! control print
360         WRITE(numout,*)
361         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
362         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
363         WRITE(numout,*) '   Namelist namicethd_pnd:'
364         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
365         WRITE(numout,*) '         Evolutive  melt pond fraction and depth                  ln_pnd_H12   = ', ln_pnd_H12
366         WRITE(numout,*) '            Melt ponds can have frozen lids                       ln_pnd_lids  = ', ln_pnd_lids
367         WRITE(numout,*) '            Allow ponds to flush thru the ice                     ln_pnd_flush = ', ln_pnd_flush
368         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
369         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
370         WRITE(numout,*) '         Prescribed melt pond fraction and depth                  ln_pnd_CST   = ', ln_pnd_CST
371         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
372         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
373         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
374      ENDIF
375      !
376      !                             !== set the choice of ice pond scheme ==!
377      ioptio = 0
378      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
379      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
380      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
381      IF( ioptio /= 1 )   &
382         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
383      !
384      SELECT CASE( nice_pnd )
385      CASE( np_pndNO )         
386         IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
387      END SELECT
388      !
389   END SUBROUTINE ice_thd_pnd_init
390   
391#else
392   !!----------------------------------------------------------------------
393   !!   Default option          Empty module          NO SI3 sea-ice model
394   !!----------------------------------------------------------------------
395#endif 
396
397   !!======================================================================
398END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.