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

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90 @ 12439

Last change on this file since 12439 was 12439, checked in by dancopsey, 4 years ago
  • Provide maximum limits to how big ponds can get. If they exceep this they leak water into the ocean.
  • Water going into melt ponds does not affect ice thickness until it leaks inot the ocean.
  • Deep snow on sea ice can cover up the ponds by forming a lid.
File size: 19.0 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    !  pond lid thickness above which the ponds disappear from the albedo calculation
42   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005    !  pond lid thickness below which the full pond area is used in the albedo calculation
43
44   REAL(wp), PARAMETER :: snow_lid_min = 0.15  ! Snow thickness above which form a lid of size pnd_lid_min on the melt ponds
45   REAL(wp), PARAMETER :: snow_lid_max = 0.25  ! Snow thickness above which form a lid of size pnd_lid_max on the melt ponds
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ice_thd_pnd
57      !!-------------------------------------------------------------------
58      !!               ***  ROUTINE ice_thd_pnd   ***
59      !!               
60      !! ** Purpose :   change melt pond fraction
61      !!               
62      !! ** Method  :   brut force
63      !!-------------------------------------------------------------------
64      !
65      SELECT CASE ( nice_pnd )
66      !
67      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
68         !
69      CASE (np_pndH12)   ;   CALL pnd_H12    !==  Holland et al 2012 melt ponds  ==!
70         !
71      END SELECT
72      !
73   END SUBROUTINE ice_thd_pnd 
74
75
76   SUBROUTINE pnd_CST 
77      !!-------------------------------------------------------------------
78      !!                ***  ROUTINE pnd_CST  ***
79      !!
80      !! ** Purpose :   Compute melt pond evolution
81      !!
82      !! ** Method  :   Melt pond fraction and thickness are prescribed
83      !!                to non-zero values when t_su = 0C
84      !!
85      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
86      !!               
87      !! ** Note   : Coupling with such melt ponds is only radiative
88      !!             Advection, ridging, rafting... are bypassed
89      !!
90      !! ** References : Bush, G.W., and Trump, D.J. (2017)
91      !!-------------------------------------------------------------------
92      INTEGER  ::   ji        ! loop indices
93      !!-------------------------------------------------------------------
94      DO ji = 1, npti
95         !
96         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
97            a_ip_frac_1d(ji) = rn_apnd
98            h_ip_1d(ji)      = rn_hpnd   
99            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
100         ELSE
101            a_ip_frac_1d(ji) = 0._wp
102            h_ip_1d(ji)      = 0._wp   
103            a_ip_1d(ji)      = 0._wp
104         ENDIF
105         !
106      END DO
107      !
108   END SUBROUTINE pnd_CST
109
110
111   SUBROUTINE pnd_H12
112      !!-------------------------------------------------------------------
113      !!                ***  ROUTINE pnd_H12  ***
114      !!
115      !! ** Purpose    : Compute melt pond evolution
116      !!
117      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds
118      !!                 and sent to ocean when surface is freezing
119      !!
120      !!                 pond growth:      Vp = Vp + dVmelt
121      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
122      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
123      !!                    with Tp = -2degC
124      !! 
125      !! ** Tunable parameters : (no real expertise yet, ideas?)
126      !!
127      !! ** Note       : Stolen from CICE for quick test of the melt pond
128      !!                 radiation and freshwater interfaces
129      !!                 Coupling can be radiative AND freshwater
130      !!                 Advection, ridging, rafting are called
131      !!
132      !! ** References : Holland, M. M. et al (J Clim 2012)
133      !!-------------------------------------------------------------------
134      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
135      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            -
136      REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_wp   ! pond aspect ratio
137      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
138      !
139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
140      REAL(wp) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness 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_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean
152      REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min
153      !
154      INTEGER  ::   ji   ! loop indices
155      !!-------------------------------------------------------------------
156      z1_rhow        = 1._wp / rhow 
157      z1_zpnd_aspect = 1._wp / zpnd_aspect
158      z1_Tp          = 1._wp / zTp 
159      z1_rhoi        = 1._wp / rhoi
160
161      ! Define time-independent field for use in refreezing
162      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
163
164      DO ji = 1, npti
165
166         IF (to_print(ji) == 10) THEN
167            write(numout,*)'icethd_pnd_start: h_ip_1d(ji), a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), ' ', a_ip_frac_1d(ji), ' ', a_ip_1d(ji)
168         END IF
169
170         !                                                        !--------------------------------!
171         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
172            !                                                     !--------------------------------!
173            !--- Remove ponds on thin ice
174            a_ip_1d(ji)      = 0._wp
175            a_ip_frac_1d(ji) = 0._wp
176            h_ip_1d(ji)      = 0._wp
177            lh_ip_1d(ji)     = 0._wp
178
179            zdh_mlt = 0._wp
180            zdh_frz = 0._wp
181            !                                                     !--------------------------------!
182         ELSE                                                     ! Case ice thickness >= rn_himin !
183            !                                                     !--------------------------------!
184            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
185            !
186            ! available meltwater for melt ponding
187            ! This is the change in ice thickness due to melt scaled up by the realive areas of the meltpond
188            ! and the area of sea ice feeding the melt ponds.
189            zdh_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) / a_ip_1d(ji)
190            !
191            !--- Pond gowth ---!
192            v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji)
193            !
194            !--- Lid shrinking. ---!
195            lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_mlt
196            !
197            !
198            !--- Pond contraction (due to refreezing) ---!
199            IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
200               t_grad = (zTp+rt0) - t_su_1d(ji)
201               
202               ! The following equation is a rearranged form of:
203               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
204               ! where: lid_thickness_start = lh_ip_1d(ji)
205               !        lid_thickness_end = lh_ip_end
206               !        omega_dt is a bunch of terms in the equation that do not change
207               ! 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
208               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
209                             
210               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
211               zdh_frz = lh_ip_end - lh_ip_1d(ji)
212
213               ! Pond shrinking
214               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
215
216               ! Lid growing
217               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
218            ELSE
219               zdh_frz = 0._wp
220            END IF
221
222            !
223            ! Make sure pond volume or lid thickness has not gone negative
224            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
225            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
226            !
227            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
228            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
229            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
230            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
231            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
232
233            ! If lid thickness is ten times greater than pond thickness then remove pond           
234            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
235                a_ip_1d(ji)      = 0._wp
236                a_ip_frac_1d(ji) = 0._wp
237                h_ip_1d(ji)      = 0._wp
238                lh_ip_1d(ji)     = 0._wp
239                v_ip_1d(ji)      = 0._wp
240            END IF
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                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji)
246                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji)
247                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
248                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
249
250                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
251                ! Scale this up to make water depth meaned over sea ice.
252                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
253
254                ! Output this water loss as a mass flux diagnostic
255                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
256
257                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
258                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
259            ENDIF
260
261            ! If pond depth exceeds half the ice thickness then reduce the pond volume
262            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
263                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
264                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji)
265                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
266                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
267                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
268
269                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
270                ! Scale this up to make water depth meaned over sea ice.
271                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
272
273                ! Output this water loss as a mass flux diagnostic
274                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
275
276                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
277                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
278            ENDIF
279
280            ! If any of the previous two IF tests has removed all the ice thickness then remove ice area.
281            IF ( h_i_1d(ji) == 0._wp ) THEN
282                a_i_1d(ji) = 0._wp
283                h_s_1d(ji) = 0._wp
284            ENDIF
285
286            ! If snow thickness exceeds snow_lid_min then form a very thin lid (so snow can go over the top)
287            IF ( h_s_1d(ji) > snow_lid_min .AND. h_s_1d(ji) < snow_lid_max ) THEN
288               weighted_lid_snow = (snow_lid_max - h_s_1d(ji))/(snow_lid_max-snow_lid_min) * pnd_lid_min +             &
289                                   (h_s_1d(ji) - snow_lid_min)/(snow_lid_max-snow_lid_min) * pnd_lid_max
290               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), weighted_lid_snow)
291            ENDIF
292            IF ( h_s_1d(ji) >= snow_lid_max ) THEN
293               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), pnd_lid_max)
294            ENDIF
295
296            !
297         ENDIF
298
299         IF (to_print(ji) == 10) THEN
300            write(numout,*)'icethd_pnd: h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji)
301            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdh_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdh_mlt
302            write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
303            write(numout,*)'icethd_pnd: lh_ip_1d(ji), zdh_mlt, zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', zdh_mlt, ' ', zdh_frz, ' ', t_su_1d(ji)
304            write(numout,*)'icethd_pnd: a_pnd_avail_1d(ji), at_i_1d(ji), wfx_pnd_1d(ji), h_i_1d(ji) = ', a_pnd_avail_1d(ji), ' ', at_i_1d(ji), ' ', wfx_pnd_1d(ji), ' ', h_i_1d(ji)
305         END IF
306
307      END DO
308      !
309   END SUBROUTINE pnd_H12
310
311
312   SUBROUTINE ice_thd_pnd_init 
313      !!-------------------------------------------------------------------
314      !!                  ***  ROUTINE ice_thd_pnd_init   ***
315      !!
316      !! ** Purpose : Physical constants and parameters linked to melt ponds
317      !!              over sea ice
318      !!
319      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond 
320      !!               parameter values called at the first timestep (nit000)
321      !!
322      !! ** input   :   Namelist namthd_pnd 
323      !!-------------------------------------------------------------------
324      INTEGER  ::   ios, ioptio   ! Local integer
325      !!
326      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
327      !!-------------------------------------------------------------------
328      !
329      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds 
330      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
331901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
332      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
333      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
334902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
335      IF(lwm) WRITE ( numoni, namthd_pnd )
336      !
337      IF(lwp) THEN                        ! control print
338         WRITE(numout,*)
339         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
340         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
341         WRITE(numout,*) '   Namelist namicethd_pnd:'
342         WRITE(numout,*) '      Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
343         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
344         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
345         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
346         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
347      ENDIF
348      !
349      !                             !== set the choice of ice pond scheme ==!
350      ioptio = 0
351                                                            nice_pnd = np_pndNO
352      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
353      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
354      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 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      END SELECT
360      !
361   END SUBROUTINE ice_thd_pnd_init
362   
363#else
364   !!----------------------------------------------------------------------
365   !!   Default option          Empty module          NO SI3 sea-ice model
366   !!----------------------------------------------------------------------
367#endif 
368
369   !!======================================================================
370END MODULE icethd_pnd 
Note: See TracBrowser for help on using the repository browser.