Changeset 12449 for NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

Timestamp:

2020-02-24T18:04:40+01:00 (4 years ago)

Author:

dancopsey

Message:

• Snow into meltponds contributes to snow to ice diagnostics affecting thickness and salinity of the ice.
• Avoid divide by zero errors by:
  • Removing meltponds from small ice areas.
  • Removing divide by pond fractions apart from lid shrinking there a minimum pond fraction is imnposed
• Ponds leaking into the ocean have a salinity flux to the ocean (also a new diagnostic).
• Add vertical fluxing of pond water.

File:

• NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90 (modified) (10 diffs)

NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

@@ -39 +39 @@
 
    REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp   ! Fraction of sea ice available for melt ponding
-   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015    !  pond lid thickness above which the ponds disappear from the albedo calculation
-   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005    !  pond lid thickness below which the full pond area is used in the albedo calculation
-
-   REAL(wp), PARAMETER :: snow_lid_min = 0.15  ! Snow thickness above which form a lid of size pnd_lid_min on the melt ponds
-   REAL(wp), PARAMETER :: snow_lid_max = 0.25  ! Snow thickness above which form a lid of size pnd_lid_max on the melt ponds
+   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation
+   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
+
+   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
 
    !! * Substitutions
@@ -138 +137 @@
       !
       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
-      REAL(wp) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness change)
+      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change)
       REAL(wp) ::   z1_Tp            ! inverse reference temperature
       REAL(wp) ::   z1_rhow          ! inverse freshwater density
@@ -151 +150 @@
       REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean
       REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min
-      !
-      INTEGER  ::   ji   ! loop indices
+      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond
+      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice
+      REAL(wp) ::   Sbr              ! Brine salinity
+      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction
+      REAL(wp) ::   perm             ! Permeability of the sea ice
+      REAL(wp) ::   drain            ! Amount of melt pond draining the sea ice per m2
+      REAL(wp) ::   za_ip            ! Temporary pond fraction
+ 
+      !
+      INTEGER  ::   ji, jk   ! loop indices
       !!-------------------------------------------------------------------
       z1_rhow        = 1._wp / rhow 
@@ -168 +175 @@
          END IF
 
-         !                                                        !--------------------------------!
-         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
-            !                                                     !--------------------------------!
-            !--- Remove ponds on thin ice
+         !                                                            !----------------------------------------------------!
+         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
+            !                                                         !----------------------------------------------------!
+            !--- Remove ponds on thin ice or tiny ice fractions
             a_ip_1d(ji)      = 0._wp
             a_ip_frac_1d(ji) = 0._wp
@@ -177 +184 @@
             lh_ip_1d(ji)     = 0._wp
 
-            zdh_mlt = 0._wp
+            zdv_mlt = 0._wp
             zdh_frz = 0._wp
             !                                                     !--------------------------------!
@@ -183 +190 @@
             !                                                     !--------------------------------!
             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
+
+            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06
+            za_ip = a_ip_1d(ji)
+            IF ( za_ip < epsi06 ) za_ip = epsi06
             !
             ! available meltwater for melt ponding
-            ! This is the change in ice thickness due to melt scaled up by the realive areas of the meltpond
-            ! and the area of sea ice feeding the melt ponds.
-            zdh_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) / a_ip_1d(ji)
+            ! This is the change in ice volume due to melt
+            zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji)
             !
             !--- Pond gowth ---!
-            v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji)
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
             !
             !--- Lid shrinking. ---!
-            lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_mlt
+            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
             !
             !
@@ -231 +241 @@
             h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
 
+            ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
+            IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
+                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji)
+                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji)
+                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            ENDIF
+
+            ! If pond depth exceeds half the ice thickness then reduce the pond volume
+            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
+                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji)
+                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
+                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
+                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            ENDIF
+
+            !-- Vertical flushing of pond water --!
+            ! 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.
+            ! This assumes the pond is sitting on top of the ice.
+            h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow
+
+            ! The depth through which water percolates is the distance under the melt pond
+            h_percolation = h_i_1d(ji) - h_ip_1d(ji)
+
+            ! Calculate the permeability of the ice (Assur 1958)
+            DO jk = 1, nlay_i
+                Sbr = - 1.2_wp                         &
+                      - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    &
+                      - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 &
+                      - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3
+                phi(jk) = sz_i_1d(ji,jk)/Sbr
+            END DO
+            perm = 3.0e-08_wp * (minval(phi))**3
+
+            ! Do the drainage using Darcy's law
+            IF ( perm > 0._wp ) THEN
+                drain = perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)
+
+                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                h_ip_1d(ji) = h_ip_1d(ji) - drain
+                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
+                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
+                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            ENDIF
+
             ! If lid thickness is ten times greater than pond thickness then remove pond            
             IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
@@ -240 +336 @@
             END IF
 
-            ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
-            IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
-                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
-                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji)
-                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
-
-                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
-                ! Scale this up to make water depth meaned over sea ice.
-                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
-
-                ! Output this water loss as a mass flux diagnostic
-                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
-
-                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
-                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
-            ENDIF
-
-            ! If pond depth exceeds half the ice thickness then reduce the pond volume
-            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
-                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
-                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji)
-                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
-                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-
-                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
-                ! Scale this up to make water depth meaned over sea ice.
-                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
-
-                ! Output this water loss as a mass flux diagnostic
-                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
-
-                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
-                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
-            ENDIF
-
-            ! If any of the previous two IF tests has removed all the ice thickness then remove ice area.
+            ! If any of the previous changes has removed all the ice thickness then remove ice area.
             IF ( h_i_1d(ji) == 0._wp ) THEN
                 a_i_1d(ji) = 0._wp
@@ -284 +342 @@
             ENDIF
 
-            ! If snow thickness exceeds snow_lid_min then form a very thin lid (so snow can go over the top)
-            IF ( h_s_1d(ji) > snow_lid_min .AND. h_s_1d(ji) < snow_lid_max ) THEN
-               weighted_lid_snow = (snow_lid_max - h_s_1d(ji))/(snow_lid_max-snow_lid_min) * pnd_lid_min +             &
-                                   (h_s_1d(ji) - snow_lid_min)/(snow_lid_max-snow_lid_min) * pnd_lid_max
-               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), weighted_lid_snow)
-            ENDIF
-            IF ( h_s_1d(ji) >= snow_lid_max ) THEN
-               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), pnd_lid_max)
-            ENDIF
-
             !
          ENDIF
@@ -299 +347 @@
          IF (to_print(ji) == 10) THEN
             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)
-            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
+            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt
             write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
-            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)
+            write(numout,*)'icethd_pnd: lh_ip_1d(ji), sfx_pnd_1d(ji), zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', sfx_pnd_1d(ji), ' ', zdh_frz, ' ', t_su_1d(ji)
             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)
+            write(numout,*)'icethd_pnd: drain, perm, h_percolation, h_gravity_head = ',drain, ' ', perm, ' ', h_percolation, ' ', h_gravity_head
+            write(numout,*)'icethd_pnd: t_i_1d(ji,1), sz_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', sz_i_1d(ji,1)
          END IF