Changeset 12484 for NEMO/branches/UKMO
- Timestamp:
- 2020-02-28T11:57:10+01:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE
- Files:
-
- 4 edited
-
ice1d.F90 (modified) (3 diffs)
-
icethd.F90 (modified) (1 diff)
-
icethd_dh.F90 (modified) (7 diffs)
-
icethd_pnd.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90
r12402 r12484 111 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_tot !: Snow accretion/ablation [m] 112 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_sum !: Ice surface ablation [m] 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_pnd !: Ice surface adlation into ponds [m] 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_pnd !: Snow surface adlation into ponds [m] 113 115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_itm !: Ice internal ablation [m] 114 116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bom !: Ice bottom ablation [m] … … 129 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_frac_1d !: 130 132 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: lh_ip_1d !: Ice pond lid thickness [m] 133 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_pnd_avail_1d !: Fraction of sea ice available for melt ponding 131 134 132 135 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 211 214 & dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new (jpij) , & 212 215 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , lh_ip_1d(jpij) , & 213 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , & 216 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , dh_i_pnd(jpij) , a_pnd_avail_1d(jpij) , & 217 & dh_s_pnd(jpij) , & 214 218 & sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d (jpij) , STAT=ierr(ii) ) 215 219 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90
r12402 r12484 220 220 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 221 221 dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 222 dh_i_pnd (1:npti) = 0._wp 223 dh_s_pnd (1:npti) = 0._wp 222 224 ! 223 225 CALL ice_thd_zdf ! --- Ice-Snow temperature --- ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90
r11715 r12484 24 24 USE lib_mpp ! MPP library 25 25 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 USE icethd_pnd, ONLY : a_pnd_avail 26 27 27 28 IMPLICIT NONE … … 98 99 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 99 100 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 101 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah_pnd ! Change in thickness due to melting into ponds 100 102 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 101 103 … … 112 114 CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile 113 115 END SELECT 116 117 ! Initialise fraction of sea ice available for melt ponding 118 DO ji = 1, npti 119 a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji) 120 END DO 114 121 115 122 ! initialize layer thicknesses and enthalpies … … 242 249 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 243 250 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 251 zdeltah_pnd(ji,jk) = zdeltah (ji,jk) * a_pnd_avail_1d(ji) 244 252 245 253 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice ! heat used to melt snow(W.m-2, >0) … … 248 256 ! updates available heat + thickness 249 257 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 258 dh_s_pnd(ji) = dh_s_pnd(ji) + zdeltah_pnd(ji,jk) ! Cumulate surface melt on areas contributing to melt ponds 250 259 zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 251 260 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) … … 343 352 344 353 zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0] 354 zdeltah_pnd(ji,jk) = - zfmdt * r1_rhoi * a_pnd_avail_1d(ji) ! Melt of layer jk [m, <0] into meltponds 345 355 346 356 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] … … 349 359 350 360 dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt 361 dh_i_pnd(ji) = dh_i_pnd(ji) + zdeltah_pnd(ji,jk) ! Cumulate surface melt on areas contributing to melt ponds 351 362 352 363 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90
r12402 r12484 37 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant pond scheme 38 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 39 40 REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp ! Fraction of sea ice available for melt ponding 41 REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 42 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 43 44 REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 39 45 40 46 !! * Substitutions … … 129 135 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 136 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 REAL(wp), PARAMETER :: freezing_t = 273.0_wp ! temperature below which refreezing occurs137 REAL(wp), PARAMETER :: max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1) 132 138 ! 133 139 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 134 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 135 REAL(wp) :: actual_mlt ! actual melt water used to make melt ponds (per m2). 136 ! Need to multiply this by ice area to work on volumes. 140 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) 137 141 REAL(wp) :: z1_Tp ! inverse reference temperature 138 142 REAL(wp) :: z1_rhow ! inverse freshwater density 139 143 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 144 REAL(wp) :: z1_rhoi ! inverse ice density 140 145 REAL(wp) :: zfac, zdum 141 146 REAL(wp) :: t_grad ! Temperature deficit for refreezing 142 147 REAL(wp) :: omega_dt ! Time independent accumulated variables used for freezing 143 148 REAL(wp) :: lh_ip_end ! Lid thickness at end of timestep (temporary variable) 144 REAL(wp) :: actual_frz ! Amount of melt pond that freezes 145 ! 146 INTEGER :: ji ! loop indices 149 REAL(wp) :: zdh_frz ! Amount of melt pond that freezes (m) 150 REAL(wp) :: v_ip_old ! Pond volume before leaking back to the ocean 151 REAL(wp) :: dh_ip_over ! Pond thickness change due to overflow or leaking 152 REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean 153 REAL(wp) :: weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min 154 REAL(wp) :: h_gravity_head ! Height above sea level of the top of the melt pond 155 REAL(wp) :: h_percolation ! Distance between the base of the melt pond and the base of the sea ice 156 REAL(wp) :: Sbr ! Brine salinity 157 REAL(wp), DIMENSION(nlay_i) :: phi ! liquid fraction 158 REAL(wp) :: perm ! Permeability of the sea ice 159 REAL(wp) :: za_ip ! Temporary pond fraction 160 REAL(wp) :: max_h_diff_ts ! Maximum meltpond depth change due to leaking or overflow (m per ts) 161 162 ! 163 INTEGER :: ji, jk ! loop indices 147 164 !!------------------------------------------------------------------- 148 165 z1_rhow = 1._wp / rhow 149 166 z1_zpnd_aspect = 1._wp / zpnd_aspect 150 167 z1_Tp = 1._wp / zTp 168 z1_rhoi = 1._wp / rhoi 169 max_h_diff_ts = max_h_diff_s * rdt_ice 151 170 152 171 ! Define time-independent field for use in refreezing … … 154 173 155 174 DO ji = 1, npti 156 ! !--------------------------------! 157 IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! 158 ! !--------------------------------! 159 !--- Remove ponds on thin ice 175 176 ! !----------------------------------------------------! 177 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 178 ! !----------------------------------------------------! 179 !--- Remove ponds on thin ice or tiny ice fractions 160 180 a_ip_1d(ji) = 0._wp 161 181 a_ip_frac_1d(ji) = 0._wp … … 163 183 lh_ip_1d(ji) = 0._wp 164 184 165 actual_mlt = 0._wp166 actual_frz = 0._wp185 zdv_mlt = 0._wp 186 zdh_frz = 0._wp 167 187 ! !--------------------------------! 168 188 ELSE ! Case ice thickness >= rn_himin ! 169 189 ! !--------------------------------! 170 190 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 171 ! 172 ! available meltwater for melt ponding [m, >0] and fraction 173 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 174 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 175 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 176 actual_mlt = zfr_mlt * zdv_mlt 191 192 ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 193 za_ip = a_ip_1d(ji) 194 IF ( za_ip < epsi06 ) za_ip = epsi06 195 ! 196 ! available meltwater for melt ponding 197 ! This is the change in ice volume due to melt 198 zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) 177 199 ! 178 200 !--- Pond gowth ---! 179 v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji)180 ! 181 !--- Lid shrinking ---!182 lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt201 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 202 ! 203 !--- Lid shrinking. ---! 204 lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 183 205 ! 184 206 ! 185 207 !--- Pond contraction (due to refreezing) ---! 186 IF ( t_su_1d(ji) < freezing_t.AND. v_ip_1d(ji) > 0._wp ) THEN187 t_grad = freezing_t- t_su_1d(ji)208 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 209 t_grad = (zTp+rt0) - t_su_1d(ji) 188 210 189 211 ! The following equation is a rearranged form of: … … 196 218 197 219 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 198 actual_frz = lh_ip_end - lh_ip_1d(ji)220 zdh_frz = lh_ip_end - lh_ip_1d(ji) 199 221 200 222 ! Pond shrinking 201 v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)223 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 202 224 203 225 ! Lid growing 204 lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz226 lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 205 227 ELSE 206 actual_frz = 0._wp228 zdh_frz = 0._wp 207 229 END IF 208 230 209 ! melt pond mass flux (<0)210 IF( zdv_mlt > 0._wp ) THEN211 zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice212 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac213 !214 ! adjust ice/snow melting flux to balance melt pond flux (>0)215 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )216 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)217 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum)218 ENDIF219 231 ! 220 232 ! Make sure pond volume or lid thickness has not gone negative … … 227 239 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 228 240 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 241 242 ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 243 IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 244 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 245 dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji) ! This will be a negative number 246 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 247 h_ip_1d(ji) = h_ip_1d(ji) + dh_ip_over 248 a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji) 249 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 250 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 251 252 ENDIF 253 254 ! If pond depth exceeds half the ice thickness then reduce the pond volume 255 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 256 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 257 dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number 258 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 259 h_ip_1d(ji) = h_ip_1d(ji) + dh_ip_over 260 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 261 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 262 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 263 264 ENDIF 265 266 !-- Vertical flushing of pond water --! 267 ! 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. 268 ! This assumes the pond is sitting on top of the ice. 269 h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 270 271 ! The depth through which water percolates is the distance under the melt pond 272 h_percolation = h_i_1d(ji) - h_ip_1d(ji) 273 274 ! Calculate the permeability of the ice (Assur 1958) 275 DO jk = 1, nlay_i 276 Sbr = - 1.2_wp & 277 - 21.8_wp * (t_i_1d(ji,jk)-rt0) & 278 - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 & 279 - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3 280 phi(jk) = sz_i_1d(ji,jk)/Sbr 281 END DO 282 perm = 3.0e-08_wp * (minval(phi))**3 283 284 ! Do the drainage using Darcy's law 285 IF ( perm > 0._wp ) THEN 286 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number 287 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative 288 289 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 290 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 291 h_ip_1d(ji) = h_ip_1d(ji) + dh_ip_over 292 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 293 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 294 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 295 296 ENDIF 229 297 230 298 ! If lid thickness is ten times greater than pond thickness then remove pond … … 236 304 v_ip_1d(ji) = 0._wp 237 305 END IF 306 307 ! If any of the previous changes has removed all the ice thickness then remove ice area. 308 IF ( h_i_1d(ji) == 0._wp ) THEN 309 a_i_1d(ji) = 0._wp 310 h_s_1d(ji) = 0._wp 311 ENDIF 312 238 313 ! 239 314 ENDIF 315 316 END DO 240 317 END DO 241 318 !
Note: See TracChangeset
for help on using the changeset viewer.
