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/dev_r13383_HPC-02_Daley_Tiling/src/ICE – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icethd_pnd.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 4 years ago

Merge in trunk up to [13550]

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