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

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icethd_pnd.F90 @ 13466

Last change on this file since 13466 was 13466, checked in by smasson, 4 years ago

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