Changeset 12423

Timestamp:

2020-02-20T15:29:25+01:00 (3 years ago)

Author:

dancopsey

Message:

Correct pond melting.

File:

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

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

@@ -129 +129 @@
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_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) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness change)
       REAL(wp) ::   z1_Tp            ! inverse reference temperature
       REAL(wp) ::   z1_rhow          ! inverse freshwater density
@@ -142 +139 @@
       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
+      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m)
       !
       INTEGER  ::   ji   ! loop indices
@@ -168 +165 @@
             lh_ip_1d(ji)     = 0._wp
 
-            actual_mlt = 0._wp
-            actual_frz = 0._wp
+            zdh_mlt = 0._wp
+            zdh_frz = 0._wp
             !                                                     !--------------------------------!
          ELSE                                                     ! Case ice thickness >= rn_himin !
@@ -175 +172 @@
             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
             !
-            ! available meltwater for melt ponding [m, >0] and fraction
-            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
+            ! 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_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * (zrmax * a_i_1d(ji)) / a_ip_1d(ji)
             !
             !--- Pond gowth ---!
-            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
+            v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji)
+            !
+            !--- Lid shrinking. ---!
+            lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_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)
+            IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
+               t_grad = (zTp+rt0) - t_su_1d(ji)
                
                ! The following equation is a rearranged form of:
@@ -201 +197 @@
                              
                lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
-               actual_frz = lh_ip_end - lh_ip_1d(ji)
+               zdh_frz = lh_ip_end - lh_ip_1d(ji)
 
                ! Pond shrinking
-               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)
+               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
 
                ! Lid growing
-               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz
+               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
             ELSE
-               actual_frz = 0._wp
+               zdh_frz = 0._wp
             END IF
 
             ! melt pond mass flux (<0)
-            IF( zdv_mlt > 0._wp ) THEN
-               zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice
+            IF( zdh_mlt > 0._wp ) THEN
+               zfac = zdh_mlt * zrmax * a_i_1d(ji) * rhow * r1_rdtice
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
@@ -246 +242 @@
          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, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt
+            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: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
-            write(numout,*)'icethd_pnd: lh_ip_1d(ji), actual_mlt, actual_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', actual_mlt, ' ', actual_frz, ' ', t_su_1d(ji)
+            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)
          END IF