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

Last change on this file since 12811 was 12811, checked in by clem, 13 months ago

debug a restartability issue. All sette tests passed now

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