Changeset 12401 for NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90

Timestamp:

2020-02-18T16:21:31+01:00 (4 years ago)

Author:

dancopsey

Message:

Add melt pond lid code.

File:

• NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90 (modified) (6 diffs)

NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90

@@ -129 +129 @@
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
       REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
+      REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs
       !
       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
+      REAL(wp) ::   actual_mlt       ! actual melt water used to make melt ponds (per m2).
+                                     ! Need to multiply this by ice area to work on volumes.
       REAL(wp) ::   z1_Tp            ! inverse reference temperature
       REAL(wp) ::   z1_rhow          ! inverse freshwater density
       REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
       REAL(wp) ::   zfac, zdum
+      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing
+      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing
+      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable)
+      REAL(wp) ::   actual_frz       ! Amount of melt pond that freezes
       !
       INTEGER  ::   ji   ! loop indices
@@ -142 +149 @@
       z1_zpnd_aspect = 1._wp / zpnd_aspect
       z1_Tp          = 1._wp / zTp 
+
+      ! Define time-independent field for use in refreezing
+      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
 
       DO ji = 1, npti
@@ -151 +161 @@
             a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp
+            lh_ip_1d(ji)     = 0._wp
+
+            actual_mlt = 0._wp
+            actual_frz = 0._wp
             !                                                     !--------------------------------!
          ELSE                                                     ! Case ice thickness >= rn_himin !
@@ -157 +171 @@
             !
             ! available meltwater for melt ponding [m, >0] and fraction
-            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
+            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow
             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
+            actual_mlt = zfr_mlt * zdv_mlt
             !
             !--- Pond gowth ---!
-            ! v_ip should never be negative, otherwise code crashes
-            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
-            !
+            v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji)
+            !
+            !--- Lid shrinking ---!
+            lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt
+            !
+            !
+            !--- Pond contraction (due to refreezing) ---!
+            IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN
+               t_grad = freezing_t - t_su_1d(ji)
+               
+               ! The following equation is a rearranged form of:
+               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
+               ! where: lid_thickness_start = lh_ip_1d(ji)
+               !        lid_thickness_end = lh_ip_end
+               !        omega_dt is a bunch of terms in the equation that do not change
+               ! 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
+               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
+                             
+               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
+               actual_frz = lh_ip_end - lh_ip_1d(ji)
+
+               ! Pond shrinking
+               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)
+
+               ! Lid growing
+               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz
+            ELSE
+               actual_frz = 0._wp
+            END IF
+
             ! melt pond mass flux (<0)
             IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+               zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
@@ -176 +218 @@
             ENDIF
             !
-            !--- Pond contraction (due to refreezing) ---!
-            v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
+            ! Make sure pond volume or lid thickness has not gone negative
+            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
+            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
             !
             ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
@@ -184 +227 @@
             a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
             h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
+
+            ! 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
+                a_ip_1d(ji)      = 0._wp
+                a_ip_frac_1d(ji) = 0._wp
+                h_ip_1d(ji)      = 0._wp
+                lh_ip_1d(ji)     = 0._wp
+                v_ip_1d(ji)      = 0._wp
+            END IF
             !
          ENDIF