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

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