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/UKMO/NEMO_4.0.1_fix_cpl/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90 @ 12497

Last change on this file since 12497 was 12497, checked in by dancopsey, 4 years ago

Fix volume and area calculations after melt ponds overflow to the ocean.
File size: 19.1 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   REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp   ! Fraction of sea ice available for melt ponding
41   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation
42   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp !  pond lid thickness below which the full pond area is used in the albedo calculation
43
44   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE ice_thd_pnd
56      !!-------------------------------------------------------------------
57      !!               ***  ROUTINE ice_thd_pnd   ***
58      !!               
59      !! ** Purpose :   change melt pond fraction
60      !!               
61      !! ** Method  :   brut force
62      !!-------------------------------------------------------------------
63      !
64      SELECT CASE ( nice_pnd )
65      !
66      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
67         !
68      CASE (np_pndH12)   ;   CALL pnd_H12    !==  Holland et al 2012 melt ponds  ==!
69         !
70      END SELECT
71      !
72   END SUBROUTINE ice_thd_pnd 
73
74
75   SUBROUTINE pnd_CST 
76      !!-------------------------------------------------------------------
77      !!                ***  ROUTINE pnd_CST  ***
78      !!
79      !! ** Purpose :   Compute melt pond evolution
80      !!
81      !! ** Method  :   Melt pond fraction and thickness are prescribed
82      !!                to non-zero values when t_su = 0C
83      !!
84      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
85      !!               
86      !! ** Note   : Coupling with such melt ponds is only radiative
87      !!             Advection, ridging, rafting... are bypassed
88      !!
89      !! ** References : Bush, G.W., and Trump, D.J. (2017)
90      !!-------------------------------------------------------------------
91      INTEGER  ::   ji        ! loop indices
92      !!-------------------------------------------------------------------
93      DO ji = 1, npti
94         !
95         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
96            a_ip_frac_1d(ji) = rn_apnd
97            h_ip_1d(ji)      = rn_hpnd   
98            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
99         ELSE
100            a_ip_frac_1d(ji) = 0._wp
101            h_ip_1d(ji)      = 0._wp   
102            a_ip_1d(ji)      = 0._wp
103         ENDIF
104         !
105      END DO
106      !
107   END SUBROUTINE pnd_CST
108
109
110   SUBROUTINE pnd_H12
111      !!-------------------------------------------------------------------
112      !!                ***  ROUTINE pnd_H12  ***
113      !!
114      !! ** Purpose    : Compute melt pond evolution
115      !!
116      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds
117      !!                 and sent to ocean when surface is freezing
118      !!
119      !!                 pond growth:      Vp = Vp + dVmelt
120      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
121      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
122      !!                    with Tp = -2degC
123      !! 
124      !! ** Tunable parameters : (no real expertise yet, ideas?)
125      !!
126      !! ** Note       : Stolen from CICE for quick test of the melt pond
127      !!                 radiation and freshwater interfaces
128      !!                 Coupling can be radiative AND freshwater
129      !!                 Advection, ridging, rafting are called
130      !!
131      !! ** References : Holland, M. M. et al (J Clim 2012)
132      !!-------------------------------------------------------------------
133      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
134      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
135      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
136      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
137      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1)
138      !
139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
140      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change)
141      REAL(wp) ::   z1_Tp            ! inverse reference temperature
142      REAL(wp) ::   z1_rhow          ! inverse freshwater density
143      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
144      REAL(wp) ::   z1_rhoi          ! inverse ice density
145      REAL(wp) ::   zfac, zdum
146      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing
147      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing
148      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable)
149      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m)
150      REAL(wp) ::   v_ip_old         ! Pond volume before leaking back to the ocean
151      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking
152      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean
153      REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min
154      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond
155      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice
156      REAL(wp) ::   Sbr              ! Brine salinity
157      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction
158      REAL(wp) ::   perm             ! Permeability of the sea ice
159      REAL(wp) ::   za_ip            ! Temporary pond fraction
160      REAL(wp) ::   max_h_diff_ts    ! Maximum meltpond depth change due to leaking or overflow (m per ts)
161     
162      !
163      INTEGER  ::   ji, jk   ! loop indices
164      !!-------------------------------------------------------------------
165      z1_rhow        = 1._wp / rhow 
166      z1_zpnd_aspect = 1._wp / zpnd_aspect
167      z1_Tp          = 1._wp / zTp 
168      z1_rhoi        = 1._wp / rhoi
169      max_h_diff_ts  = max_h_diff_s * rdt_ice
170
171      ! Define time-independent field for use in refreezing
172      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
173
174      DO ji = 1, npti
175
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            a_ip_frac_1d(ji) = 0._wp
182            h_ip_1d(ji)      = 0._wp
183            lh_ip_1d(ji)     = 0._wp
184
185            zdv_mlt = 0._wp
186            zdh_frz = 0._wp
187            !                                                     !--------------------------------!
188         ELSE                                                     ! Case ice thickness >= rn_himin !
189            !                                                     !--------------------------------!
190            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
191
192            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06
193            za_ip = a_ip_1d(ji)
194            IF ( za_ip < epsi06 ) za_ip = epsi06
195            !
196            ! available meltwater for melt ponding
197            ! This is the change in ice volume due to melt
198            zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji)
199            !
200            !--- Pond gowth ---!
201            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
202            !
203            !--- Lid shrinking. ---!
204            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
205            !
206            !
207            !--- Pond contraction (due to refreezing) ---!
208            IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
209               t_grad = (zTp+rt0) - t_su_1d(ji)
210               
211               ! The following equation is a rearranged form of:
212               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
213               ! where: lid_thickness_start = lh_ip_1d(ji)
214               !        lid_thickness_end = lh_ip_end
215               !        omega_dt is a bunch of terms in the equation that do not change
216               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same
217               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
218                             
219               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
220               zdh_frz = lh_ip_end - lh_ip_1d(ji)
221
222               ! Pond shrinking
223               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
224
225               ! Lid growing
226               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
227            ELSE
228               zdh_frz = 0._wp
229            END IF
230
231            !
232            ! Make sure pond volume or lid thickness has not gone negative
233            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
234            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
235            !
236            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
237            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
238            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
239            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
240            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
241
242            ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
243            IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
244                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
245                dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number
246                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
247                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
248                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
249                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
250                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
251
252            ENDIF
253
254            ! If pond depth exceeds half the ice thickness then reduce the pond volume
255            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
256                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
257                dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number
258                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
259                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
260                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
261                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
262                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
263
264            ENDIF
265
266            !-- Vertical flushing of pond water --!
267            ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth.
268            ! This assumes the pond is sitting on top of the ice.
269            h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow
270
271            ! The depth through which water percolates is the distance under the melt pond
272            h_percolation = h_i_1d(ji) - h_ip_1d(ji)
273
274            ! Calculate the permeability of the ice (Assur 1958)
275            DO jk = 1, nlay_i
276                Sbr = - 1.2_wp                         &
277                      - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    &
278                      - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 &
279                      - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3
280                phi(jk) = sz_i_1d(ji,jk)/Sbr
281            END DO
282            perm = 3.0e-08_wp * (minval(phi))**3
283
284            ! Do the drainage using Darcy's law
285            IF ( perm > 0._wp ) THEN
286                dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number
287                dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative
288
289                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
290                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
291                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
292                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
293                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
294                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
295
296            ENDIF
297
298            ! If lid thickness is ten times greater than pond thickness then remove pond           
299            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
300                a_ip_1d(ji)      = 0._wp
301                a_ip_frac_1d(ji) = 0._wp
302                h_ip_1d(ji)      = 0._wp
303                lh_ip_1d(ji)     = 0._wp
304                v_ip_1d(ji)      = 0._wp
305            END IF
306
307            ! If any of the previous changes has removed all the ice thickness then remove ice area.
308            IF ( h_i_1d(ji) == 0._wp ) THEN
309                a_i_1d(ji) = 0._wp
310                h_s_1d(ji) = 0._wp
311            ENDIF
312
313            !
314         ENDIF
315
316      END DO
317      !
318   END SUBROUTINE pnd_H12
319
320
321   SUBROUTINE ice_thd_pnd_init 
322      !!-------------------------------------------------------------------
323      !!                  ***  ROUTINE ice_thd_pnd_init   ***
324      !!
325      !! ** Purpose : Physical constants and parameters linked to melt ponds
326      !!              over sea ice
327      !!
328      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
329      !!               parameter values called at the first timestep (nit000)
330      !!
331      !! ** input   :   Namelist namthd_pnd 
332      !!-------------------------------------------------------------------
333      INTEGER  ::   ios, ioptio   ! Local integer
334      !!
335      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
336      !!-------------------------------------------------------------------
337      !
338      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
339      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
340901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist' )
341      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
342      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
343902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
344      IF(lwm) WRITE ( numoni, namthd_pnd )
345      !
346      IF(lwp) THEN                        ! control print
347         WRITE(numout,*)
348         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
349         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
350         WRITE(numout,*) '   Namelist namicethd_pnd:'
351         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd
352         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
353         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
354         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
355         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
356         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
357      ENDIF
358      !
359      !                             !== set the choice of ice pond scheme ==!
360      ioptio = 0
361      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF
362      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
363      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
364      IF( ioptio /= 1 )   &
365         & 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)' )
366      !
367      SELECT CASE( nice_pnd )
368      CASE( np_pndNO )         
369         IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
370      END SELECT
371      !
372   END SUBROUTINE ice_thd_pnd_init
373   
374#else
375   !!----------------------------------------------------------------------
376   !!   Default option          Empty module          NO SI3 sea-ice model
377   !!----------------------------------------------------------------------
378#endif 
379
380   !!======================================================================
381END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.