Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/field_def_nemo-ice.xml (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/field_def_nemo-ice.xml (revision 12636)
@@ -253,4 +253,6 @@
+
+
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/namelist_ice_ref (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/namelist_ice_ref (revision 12636)
@@ -177,4 +177,10 @@
ln_pnd = .false. ! activate melt ponds or not
ln_pnd_H12 = .false. ! activate evolutive melt ponds (from Holland et al 2012)
+ rn_pnd_min = 0.15 ! Minimum ice fraction that contributes to melt ponds
+ rn_pnd_max = 0.7 ! Maximum ice fraction that contributes to melt ponds
+ ln_pnd_totfrac = .false. ! Use total ice fraction instead of category ice fraction in melt pond contribution fracton
+ ln_pnd_overflow = .false. ! Allow ponds to overflow and have vertical flushing
+ ln_pnd_lids = .false. ! Melt ponds can have frozen lids
+ ln_use_pndmass = .true. ! Use melt pond mass flux diagnostic, passing it to the ocean for emp
ln_pnd_CST = .false. ! activate constant melt ponds
rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/doc/namelists/namthd_pnd
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/doc/namelists/namthd_pnd (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/doc/namelists/namthd_pnd (revision 12636)
@@ -4,4 +4,10 @@
ln_pnd = .false. ! activate melt ponds or not
ln_pnd_H12 = .false. ! activate evolutive melt ponds (from Holland et al 2012)
+ rn_pnd_min = 0.15 ! Minimum ice fraction that contributes to melt ponds
+ rn_pnd_max = 0.7 ! Maximum ice fraction that contributes to melt ponds
+ ln_pnd_totfrac = .false. ! Use total ice fraction instead of category ice fraction in melt pond contribution fracton
+ ln_pnd_overflow = .false. ! Allow ponds to overflow and have vertical flushing
+ ln_pnd_lids = .false. ! Melt ponds can have frozen lids
+ ln_use_pndmass = .true. ! Use melt pond mass flux diagnostic, passing it to the ocean for emp
ln_pnd_CST = .false. ! activate constant melt ponds
rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90 (revision 12636)
@@ -70,4 +70,5 @@
!! a_ip | - | Ice pond concentration | |
!! v_ip | - | Ice pond volume per unit area| m |
+ !! lh_ip ! lh_ip_1d ! Ice pond lid thickness ! m !
!! |
!!-------------|-------------|---------------------------------|-------|
@@ -195,4 +196,11 @@
REAL(wp), PUBLIC :: rn_hpnd !: prescribed pond depth (0 0 ) THEN
+ CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip)
+ ELSE
+ lh_ip(:,:,:) = 0._wp
+ ENDIF
! fields needed for Met Office (Jules) coupling
IF( ln_cpl ) THEN
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90 (revision 12636)
@@ -132,5 +132,5 @@
! --- cloud-sky and overcast-sky ice albedos --- !
- CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
+ CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os )
! albedo depends on cloud fraction because of non-linear spectral effects
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90 (revision 12636)
@@ -355,4 +355,6 @@
CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+ CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+ CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) )
!
CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice )
@@ -461,4 +463,6 @@
CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+ CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+ CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) )
!
CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90 (revision 12636)
@@ -38,4 +38,6 @@
INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012)
+ REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
+
!! * Substitutions
# include "vectopt_loop_substitute.h90"
@@ -89,8 +91,10 @@
IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
a_ip_frac_1d(ji) = rn_apnd
+ a_ip_eff_1d(ji) = rn_apnd
h_ip_1d(ji) = rn_hpnd
a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
ELSE
a_ip_frac_1d(ji) = 0._wp
+ a_ip_eff_1d(ji) = 0._wp
h_ip_1d(ji) = 0._wp
a_ip_1d(ji) = 0._wp
@@ -129,43 +133,81 @@
REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio
REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature
- !
+ REAL(wp), PARAMETER :: max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1)
+ REAL(wp), PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation
+ REAL(wp), PARAMETER :: pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation
+ !
+ REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean)
REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding
- REAL(wp) :: zdv_mlt ! available meltwater for melt ponding
+ 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
REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio
+ REAL(wp) :: z1_rhoi ! inverse ice density
REAL(wp) :: zfac, zdum
- !
- INTEGER :: ji ! loop indices
+ 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) :: zdh_frz ! Amount of melt pond that freezes (m)
+ REAL(wp) :: v_ip_old ! Pond volume before leaking back to the ocean
+ REAL(wp) :: dh_ip_over ! Pond thickness change due to overflow or leaking
+ REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean
+ 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) :: za_ip ! Temporary pond fraction
+ REAL(wp) :: max_h_diff_ts ! Maximum meltpond depth change due to leaking or overflow (m per ts)
+ REAL(wp) :: lfrac_pnd ! The fraction of the meltpond exposed (not inder a frozen lid)
+
+ !
+ INTEGER :: ji, jk ! loop indices
!!-------------------------------------------------------------------
z1_rhow = 1._wp / rhow
z1_zpnd_aspect = 1._wp / zpnd_aspect
z1_Tp = 1._wp / zTp
+ z1_rhoi = 1._wp / rhoi
+ max_h_diff_ts = max_h_diff_s * rdt_ice
+
+ ! Define time-independent field for use in refreezing
+ omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
DO ji = 1, npti
- ! !--------------------------------!
- 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
h_ip_1d(ji) = 0._wp
+ lh_ip_1d(ji) = 0._wp
! !--------------------------------!
ELSE ! Case ice thickness >= rn_himin !
! !--------------------------------!
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 * a_i_1d(ji)
- zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc
- !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper
+
+ ! 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 [m, >0] and fraction
+ tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow
+ IF ( ln_pnd_totfrac ) THEN
+ zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction
+ ELSE
+ zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji) ! Use category ice fraction
+ ENDIF
+ zdv_mlt = zfr_mlt * tot_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) + zdv_mlt
+ !
+ !--- Lid shrinking. ---!
+ IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
!
! melt pond mass flux (<0)
- IF( zdv_mlt > 0._wp ) THEN
- zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+ IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN
+ zfac = zdv_mlt * rhow * r1_rdtice
wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
!
@@ -177,5 +219,39 @@
!
!--- 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 )
+ IF ( ln_pnd_lids ) THEN
+
+ ! Code to use if using melt pond lids
+ 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:
+ ! 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 mathematically easier
+ ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).
+
+ lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
+ zdh_frz = lh_ip_end - lh_ip_1d(ji)
+
+ ! Pond shrinking
+ v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
+
+ ! Lid growing
+ IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
+ ELSE
+ zdh_frz = 0._wp
+ END IF
+
+ ELSE
+ ! Code to use if not using melt pond lids
+ v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
+ ENDIF
+
+ !
+ ! 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,6 +260,96 @@
a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji)
+
+ !--- Pond overflow and vertical flushing ---!
+ IF ( ln_pnd_overflow ) THEN
+
+ ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+ IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN
+ v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use
+ dh_ip_over = zpnd_aspect * zfr_mlt - h_ip_1d(ji) ! This will be a negative number
+ dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit
+ h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+ 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)
+ 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
+ dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number
+ dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit
+ h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+ 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)
+ 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
+ dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number
+ dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative
+
+ v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use
+ dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit
+ h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+ 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)
+ ENDIF
+ ENDIF
+
+ ! If lid thickness is ten times greater than pond thickness then remove pond
+ IF ( ln_pnd_lids ) THEN
+ 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
+ ENDIF
+ ENDIF
+
+ ! 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
+ h_s_1d(ji) = 0._wp
+ ENDIF
+
!
ENDIF
+
+ ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid)
+ IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN ! Pond lid is very thin. Expose all the pond.
+ lfrac_pnd = 1.0
+ ELSE
+ IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN ! Pond lid is very thick. Cover all the pond up with ice and nosw.
+ lfrac_pnd = 0.0
+ ELSE ! Pond lid is of intermediate thickness. Expose part of the melt pond.
+ lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min)
+ ENDIF
+ ENDIF
+
+ ! Calculate the melt pond effective area
+ a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd
+
END DO
!
@@ -205,5 +371,7 @@
INTEGER :: ios, ioptio ! Local integer
!!
- NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+ NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb, &
+ rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac, &
+ ln_use_pndmass
!!-------------------------------------------------------------------
!
@@ -223,4 +391,10 @@
WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd
WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
+ WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_pnd_min = ', rn_pnd_min
+ WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_pnd_max = ', rn_pnd_max
+ WRITE(numout,*) ' Use total ice fraction instead of category ice fraction ln_pnd_totfrac = ',ln_pnd_totfrac
+ WRITE(numout,*) ' Allow ponds to overflow and have vertical flushing ln_pnd_overflow = ',ln_pnd_overflow
+ WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ',ln_pnd_lids
+ WRITE(numout,*) ' Use melt pond mass flux diagnostic, passing to ocean ln_use_pndmass = ',ln_use_pndmass
WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST
WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90 (revision 12636)
@@ -185,5 +185,5 @@
! Snow/ice albedo (only if sent to coupler, useless in forced mode)
!------------------------------------------------------------------
- CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+ CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
!
alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90 (revision 12636)
@@ -533,4 +533,5 @@
a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+ lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj)
!
END DO
@@ -555,5 +556,5 @@
- SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+ SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_zapneg ***
@@ -570,4 +571,5 @@
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
+ REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
@@ -637,9 +639,10 @@
WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
! but it does not change conservation, so keep it this way is ok
+ WHERE( plh_ip (:,:,:) < 0._wp ) plh_ip (:,:,:) = 0._wp
!
END SUBROUTINE ice_var_zapneg
- SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+ SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_roundoff ***
@@ -654,4 +657,5 @@
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
+ REAL(wp), DIMENSION(:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
@@ -668,4 +672,5 @@
WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
+ WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 ) plh_ip(1:npti,:) = 0._wp ! lh_ip must be >= 0
ENDIF
!
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90 (revision 12636)
@@ -149,4 +149,5 @@
! --- category-dependent fields --- !
+ IF( iom_use('icelhpnd_cat') ) CALL iom_put( 'icelhpnd_cat', lh_ip * zmsk00l ) ! melt pond lid thickness for categories
IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0%
IF( iom_use('iceconc_cat' ) ) CALL iom_put( 'iceconc_cat' , a_i * zmsk00l ) ! area for categories
@@ -165,4 +166,7 @@
IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac for categories
IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
+
+ IF( iom_use('icevol_cat' ) ) CALL iom_put( "icevol_cat" , v_i * zmsk00l ) ! volume for categories
+ IF( iom_use('iceaepnd_cat') ) CALL iom_put( 'iceaepnd_cat', a_ip_eff * zmsk00l ) ! melt pond effective frac for categories
!------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/BDY/bdyice.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/BDY/bdyice.F90 (revision 12636)
@@ -98,4 +98,5 @@
& , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1. &
& , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+ CALL lbc_lnk_multi( 'bdyice', lh_ip,'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
! exchange 4d arrays : third dimension = 1 and then third dimension = jpk
CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
@@ -163,4 +164,5 @@
a_ip(ji,jj, jl) = ( a_ip(ji,jj, jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond concentration
h_ip(ji,jj, jl) = ( h_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond depth
+ lh_ip(ji,jj, jl) = ( lh_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond lid thickness
!
sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
@@ -170,4 +172,5 @@
a_ip(ji,jj,jl) = 0._wp
h_ip(ji,jj,jl) = 0._wp
+ lh_ip(ji,jj,jl) = 0._wp
ENDIF
!
@@ -231,4 +234,5 @@
a_ip(ji,jj, jl) = a_ip(ib,jb, jl)
h_ip(ji,jj, jl) = h_ip(ib,jb, jl)
+ lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl)
!
sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
@@ -280,4 +284,5 @@
a_ip(ji,jj, jl) = 0._wp
v_ip(ji,jj, jl) = 0._wp
+ lh_ip(ji,jj, jl) = 0._wp
t_su(ji,jj, jl) = rt0
t_s (ji,jj,:,jl) = rt0
Index: /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/SBC/sbccpl.F90 (revision 12635)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/SBC/sbccpl.F90 (revision 12636)
@@ -2047,9 +2047,11 @@
INTEGER, INTENT(in) :: kt
!
+ !
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: isec, info ! local integer
REAL(wp) :: zumax, zvmax
REAL(wp), DIMENSION(jpi,jpj) :: zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
- REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztmp3, ztmp4
+ REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztmp3, ztmp4
+
!!----------------------------------------------------------------------
!
@@ -2256,6 +2258,8 @@
SELECT CASE( sn_snd_mpnd%clcat )
CASE( 'yes' )
- ztmp3(:,:,1:jpl) = a_ip(:,:,1:jpl)
- ztmp4(:,:,1:jpl) = v_ip(:,:,1:jpl)
+
+ ztmp3(:,:,1:jpl) = a_ip_eff(:,:,1:jpl)
+ ztmp4(:,:,1:jpl) = h_ip(:,:,1:jpl)
+
CASE( 'no' )
ztmp3(:,:,:) = 0.0