Changeset 8906 for branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90

Timestamp:

2017-12-05T18:24:20+01:00 (6 years ago)

Author:

clem

Message:

make melt ponds from Holland2012 and associated freshwater flux conservative

File:

• branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90 (modified) (5 diffs)

branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90

@@ -125 +125 @@
       !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
       !!                    with Tp = -2degC
-      !!
+      !!  
       !! ** Tunable parameters : (no real expertise yet, ideas?)
       !! 
@@ -136 +136 @@
       !!    
       !!-------------------------------------------------------------------
-      REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
-      REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
-      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
-      REAL(wp), PARAMETER ::   zTp = -2._wp ! 
-
-      REAL(wp) ::   zfr_mlt            ! fraction of available meltwater retained for melt ponding
+      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
+      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
+      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
+      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
+
+      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
-      REAL(wp) ::   z1_Tp        ! reference temperature
-      REAL(wp) ::   z1_rhofw          ! inverse freshwater density
-      REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
-      REAL(wp) ::   zvpnd            ! dummy pond volume
+      REAL(wp) ::   z1_Tp            ! inverse reference temperature
+      REAL(wp) ::   z1_rhofw         ! inverse freshwater density
+      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
       REAL(wp) ::   zfac, zdum
 
-      INTEGER  ::   ji        ! loop indices
+      INTEGER  ::   ji   ! loop indices
       !!-------------------------------------------------------------------
       z1_rhofw       = 1. / rhofw 
@@ -159 +158 @@
          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
             !                                                     !--------------------------------!
-            !--- Remove ponds on thin ice and send water into the ocean
-            IF( ln_pnd_fwb ) THEN
-               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice
-            ENDIF
+            !--- Remove ponds on thin ice
             a_ip_1d(ji)      = 0._wp
             a_ip_frac_1d(ji) = 0._wp
@@ -170 +166 @@
             !                                                     !--------------------------------!
             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
-            zvpnd       = v_ip_1d(ji)
 
             ! available meltwater for melt ponding [m, >0] and fraction
             zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji)
-            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  !!clem CICE doc
-!!!            zfr_mlt = zrmin + zrmax * a_i_1d(ji)  !!clem Holland paper 
-
-            !--- Pond gowth
+            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
+            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
+
+            !--- Pond gowth ---!
             ! v_ip should never be negative, otherwise code crashes
             ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather)
             v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
 
+            ! melt pond mass flux (<0)
             IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN
-               ! melt ponds mass flux (<0)
                zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
 
-               ! adjust ice/snow melting (>0) to take into account ponding
+               ! adjust ice/snow melting flux to balance melt pond flux (>0)
                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
@@ -193 +188 @@
             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 )
-
-            ! --- Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
-            !        h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
+            !--- 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 )
+
+            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
+            !    h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
             a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
-            !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
             a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
             h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
-
-            !! clem: there is no clever way to conserve mass here...
-            
-!            IF( ln_pnd_fwb ) THEN
-!               ! melt ponds mass flux (>0 when ponds shrink)
-!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice
-!               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!!
-!
-!               ! adjust ice/snow
-!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )!
-!
-!               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice
-!               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice
-!               
-!               !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac )
-!               !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )!
-!
-!            ENDIF
 
          ENDIF