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

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

clean useless variables

  • Property svn:keywords set to Id
File size: 19.6 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            ! clem: problem with conservation or not ?
188            !                                                         !--------------------------------!
189         ELSE                                                         ! Case ice thickness >= rn_himin !
190            !                                                         !--------------------------------!
191            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
192            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
193            !
194            !------------------!
195            ! case ice melting !
196            !------------------!
197            !
198            !--- available meltwater for melt ponding ---!
199            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
200            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
201            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?
202            !
203            !--- overflow ---!
204            ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
205            !    a_ip_max = zfr_mlt * a_i
206            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
207            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
208            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
209
210            ! If pond depth exceeds half the ice thickness then reduce the pond volume
211            !    h_ip_max = 0.5 * h_i
212            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
213            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
214            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
215           
216            !--- Pond growing ---!
217            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
218            !
219            !--- Lid melting ---!
220            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
221            !
222            !--- mass flux ---!
223            IF( zdv_mlt > 0._wp ) THEN
224               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
225               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
226               !
227               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
228               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
229               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
230            ENDIF
231
232            !-------------------!
233            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
234            !-------------------!
235            !
236            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
237            !
238            !--- Pond contraction (due to refreezing) ---!
239            IF( ln_pnd_lids ) THEN
240               !
241               !--- Lid growing and subsequent pond shrinking ---!
242               zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
243                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
244               
245               ! Lid growing
246               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
247               
248               ! Pond shrinking
249               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
250
251            ELSE
252               ! Pond shrinking
253               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
254            ENDIF
255            !
256            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
257            ! v_ip     = h_ip * a_ip
258            ! 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)
259            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
260            h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
261
262            !---------------!           
263            ! Pond flushing !
264            !---------------!
265            IF( ln_pnd_flush ) THEN
266               ! height of top of the pond above sea-level
267               zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
268
269               ! Calculate the permeability of the ice (Assur 1958)
270               DO jk = 1, nlay_i
271                  zsbr = - 1.2_wp                                  &
272                     &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
273                     &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
274                     &   - 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)
275                  ztmp(jk) = sz_i_1d(ji,jk) / zsbr
276               END DO
277               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
278
279               ! Do the drainage using Darcy's law
280               zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
281               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
282               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
283               
284               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
285               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
286               h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
287
288            ENDIF
289
290            !--- Corrections and lid thickness ---!
291            IF( ln_pnd_lids ) THEN
292               !--- retrieve lid thickness from volume ---!
293               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
294               ELSE                              ;   h_il_1d(ji) = 0._wp
295               ENDIF
296               !--- remove ponds if lids are much larger than ponds ---!
297               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
298                  a_ip_1d(ji)      = 0._wp
299                  h_ip_1d(ji)      = 0._wp
300                  h_il_1d(ji)      = 0._wp
301               ENDIF
302            ENDIF
303            !
304         ENDIF
305         
306      END DO
307      !
308   END SUBROUTINE pnd_H12
309
310
311   SUBROUTINE ice_thd_pnd_init 
312      !!-------------------------------------------------------------------
313      !!                  ***  ROUTINE ice_thd_pnd_init   ***
314      !!
315      !! ** Purpose : Physical constants and parameters linked to melt ponds
316      !!              over sea ice
317      !!
318      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
319      !!               parameter values called at the first timestep (nit000)
320      !!
321      !! ** input   :   Namelist namthd_pnd 
322      !!-------------------------------------------------------------------
323      INTEGER  ::   ios, ioptio   ! Local integer
324      !!
325      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, &
326         &                          ln_pnd_CST, rn_apnd, rn_hpnd, &
327         &                          ln_pnd_alb
328      !!-------------------------------------------------------------------
329      !
330      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
331      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
332901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
333      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
334      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
335902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
336      IF(lwm) WRITE ( numoni, namthd_pnd )
337      !
338      IF(lwp) THEN                        ! control print
339         WRITE(numout,*)
340         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
341         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
342         WRITE(numout,*) '   Namelist namicethd_pnd:'
343         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
344         WRITE(numout,*) '         Evolutive  melt pond fraction and depth                  ln_pnd_H12   = ', ln_pnd_H12
345         WRITE(numout,*) '            Melt ponds can have frozen lids                       ln_pnd_lids  = ', ln_pnd_lids
346         WRITE(numout,*) '            Allow ponds to flush thru the ice                     ln_pnd_flush = ', ln_pnd_flush
347         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
348         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max
349         WRITE(numout,*) '         Prescribed melt pond fraction and depth                  ln_pnd_CST   = ', ln_pnd_CST
350         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd
351         WRITE(numout,*) '            Prescribed pond depth                                 rn_hpnd      = ', rn_hpnd
352         WRITE(numout,*) '         Melt ponds affect albedo or not                          ln_pnd_alb   = ', ln_pnd_alb
353      ENDIF
354      !
355      !                             !== set the choice of ice pond scheme ==!
356      ioptio = 0
357      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
358      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
359      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
360      IF( ioptio /= 1 )   &
361         & 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)' )
362      !
363      SELECT CASE( nice_pnd )
364      CASE( np_pndNO )         
365         IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
366      END SELECT
367      !
368   END SUBROUTINE ice_thd_pnd_init
369   
370#else
371   !!----------------------------------------------------------------------
372   !!   Default option          Empty module          NO SI3 sea-ice model
373   !!----------------------------------------------------------------------
374#endif 
375
376   !!======================================================================
377END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.