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

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

fix previous revision. Make sure a_ip and a_ip_eff are correct

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