- Timestamp:
- 2020-11-28T10:38:45+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/SI3-05_MeltPonds_topo
- Files:
-
- 12 edited
-
cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml (modified) (3 diffs)
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (2 diffs)
-
cfgs/SHARED/namelist_ice_ref (modified) (2 diffs)
-
src/ICE/ice.F90 (modified) (1 diff)
-
src/ICE/icedyn_adv_pra.F90 (modified) (9 diffs)
-
src/ICE/icedyn_adv_umx.F90 (modified) (3 diffs)
-
src/ICE/icedyn_rdgrft.F90 (modified) (3 diffs)
-
src/ICE/iceistate.F90 (modified) (1 diff)
-
src/ICE/iceitd.F90 (modified) (2 diffs)
-
src/ICE/icethd_pnd.F90 (modified) (28 diffs)
-
src/ICE/icevar.F90 (modified) (2 diffs)
-
src/ICE/icewri.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml
r9896 r13902 28 28 <field field_ref="iceapnd" name="siapnd" /> 29 29 <field field_ref="icevpnd" name="sivpnd" /> 30 <field field_ref="icevlid" name="sivlid" /> 30 31 <field field_ref="sst_m" name="sst_m" /> 31 32 <field field_ref="sss_m" name="sss_m" /> … … 36 37 <field field_ref="icetbot" name="sitbot" /> 37 38 <field field_ref="icetsni" name="sitsni" /> 39 40 <!-- ponds --> 41 <field field_ref="dvpn_mlt" name="dvpn_mlt" /> 42 <field field_ref="dvpn_lid" name="dvpn_lid" /> 43 <field field_ref="dvpn_rnf" name="dvpn_rnf" /> 44 <field field_ref="dvpn_drn" name="dvpn_drn" /> 38 45 39 46 <!-- momentum --> … … 85 92 <field field_ref="icesalt_cat" name="sisalcat"/> 86 93 <field field_ref="icetemp_cat" name="sitemcat"/> 94 <field field_ref="iceapnd_cat" name="siapncat"/> 95 <field field_ref="icevpnd_cat" name="sivpncat"/> 87 96 88 97 </file> -
NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/field_def_nemo-ice.xml
r13648 r13902 51 51 <field id="icehlid" long_name="melt pond lid depth" standard_name="sea_ice_meltpondlid_depth" unit="m" /> 52 52 <field id="icevlid" long_name="melt pond lid volume" standard_name="sea_ice_meltpondlid_volume" unit="m" /> 53 54 <field id="dvpn_mlt" long_name="pond volume tendency due to surface melt" standard_name="sea_ice_pondvolume_tendency_melt" unit="cm/d" /> 55 <field id="dvpn_lid" long_name="pond volume tendency due to exchanges with lid" standard_name="sea_ice_pondvolume_tendency_lids" unit="cm/d" /> 56 <field id="dvpn_rnf" long_name="pond volume tendency due to runoff" standard_name="sea_ice_pondvolume_tendency_runoff" unit="cm/d" /> 57 <field id="dvpn_drn" long_name="pond volume tendency due to drainage" standard_name="sea_ice_pondvolume_tendency_drainage" unit="cm/d" /> 53 58 54 59 <!-- heat --> … … 295 300 <field id="snwtemp_cat" long_name="Snow temperature per category" unit="degC" detect_missing_value="true" /> 296 301 <field id="icettop_cat" long_name="Ice/snow surface temperature per category" unit="degC" detect_missing_value="true" /> 297 <field id="iceapnd_cat" long_name="Ice melt pond concentration per category" unit="" /> 302 <field id="icevpnd_cat" long_name="Ice melt pond volume per grid area per category" unit="m" /> 303 <field id="iceapnd_cat" long_name="Ice melt pond grid fraction" unit="" /> 298 304 <field id="icehpnd_cat" long_name="Ice melt pond thickness per category" unit="m" detect_missing_value="true" /> 299 305 <field id="icehlid_cat" long_name="Ice melt pond lid thickness per category" unit="m" detect_missing_value="true" /> 300 <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"unit="" />306 <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category" unit="" /> 301 307 <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category" unit="" /> 302 308 <field id="icemask_cat" long_name="Fraction of time step with sea ice (per category)" unit="" /> -
NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/namelist_ice_ref
r13850 r13902 195 195 !------------------------------------------------------------------------------ 196 196 ln_pnd = .true. ! activate melt ponds or not 197 ln_pnd_TOPO = .false. ! topographic melt ponds 198 ln_pnd_LEV = .true. ! level ice melt ponds 199 rn_apnd_min = 0.15 ! minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 200 rn_apnd_max = 0.85 ! maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 197 ln_pnd_TOPO = .true. ! topographic melt ponds 198 ln_pnd_LEV = .false. ! level ice melt ponds 199 rn_apnd_min = 0.15 ! minimum meltwater fraction contributing to pond growth (TOPO and LEV) 200 rn_apnd_max = 0.85 ! maximum meltwater fraction contributing to pond growth (TOPO and LEV) 201 rn_pnd_flush= 0.1 ! pond flushing efficiency (tuning parameter) (LEV) 201 202 ln_pnd_CST = .false. ! constant melt ponds 202 203 rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC 203 204 rn_hpnd = 0.05 ! prescribed pond depth, at Tsu=0 degC 204 ln_pnd_lids = . true.! frozen lids on top of the ponds (only for ln_pnd_LEV)205 ln_pnd_lids = .false. ! frozen lids on top of the ponds (only for ln_pnd_LEV) 205 206 ln_pnd_alb = .true. ! effect of melt ponds on ice albedo 206 207 / … … 227 228 rn_tms_ini_n = 270. ! initial snw temperature (K), North 228 229 rn_tms_ini_s = 270. ! " " South 229 rn_apd_ini_n = 0. 2! initial pond fraction (-), North230 rn_apd_ini_n = 0.0 ! initial pond fraction (-), North 230 231 rn_apd_ini_s = 0.2 ! " " South 231 rn_hpd_ini_n = 0.0 5! initial pond depth (m), North232 rn_hpd_ini_n = 0.00 ! initial pond depth (m), North 232 233 rn_hpd_ini_s = 0.05 ! " " South 233 234 rn_hld_ini_n = 0.0 ! initial pond lid depth (m), North -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/ice.F90
r13850 r13902 210 210 LOGICAL , PUBLIC :: ln_pnd_TOPO !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 211 211 LOGICAL , PUBLIC :: ln_pnd_LEV !: Simple melt pond scheme 212 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum ice fraction that contributes to melt ponds 213 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum ice fraction that contributes to melt ponds 212 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum fraction of melt water contributing to ponds 213 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum fraction of melt water contributing to ponds 214 REAL(wp), PUBLIC :: rn_pnd_flush !: Pond flushing efficiency (tuning parameter) 214 215 LOGICAL , PUBLIC :: ln_pnd_CST !: Melt ponds scheme with constant fraction and depth 215 216 REAL(wp), PUBLIC :: rn_apnd !: prescribed pond fraction (0<rn_apnd<1) -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icedyn_adv_pra.F90
r13634 r13902 178 178 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 179 179 END DO 180 IF ( ln_pnd_LEV ) THEN180 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 181 181 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 182 182 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume … … 214 214 END DO 215 215 ! 216 IF ( ln_pnd_LEV ) THEN216 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 217 217 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 218 218 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) … … 249 249 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 250 250 END DO 251 IF ( ln_pnd_LEV ) THEN251 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 252 252 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 253 253 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) … … 278 278 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy 279 279 & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp ) 280 IF ( ln_pnd_LEV ) THEN280 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 281 281 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction 282 282 & , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp & … … 302 302 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 303 303 END DO 304 IF ( ln_pnd_LEV ) THEN304 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 305 305 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 306 306 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) … … 780 780 ! ! -- check h_ip -- ! 781 781 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 782 IF( ln_pnd_LEV.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN782 IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 783 783 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 784 784 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 1027 1027 END DO 1028 1028 ! 1029 IF( ln_pnd_LEV ) THEN ! melt pond fraction1029 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction 1030 1030 IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 1031 1031 CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap ) … … 1069 1069 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 1070 1070 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content 1071 IF( ln_pnd_LEV ) THEN1071 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 1072 1072 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 1073 1073 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume … … 1137 1137 END DO 1138 1138 ! 1139 IF( ln_pnd_LEV ) THEN ! melt pond fraction1139 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction 1140 1140 CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap ) 1141 1141 CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap ) -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icedyn_adv_umx.F90
r13617 r13902 341 341 ! 342 342 !== melt ponds ==! 343 IF ( ln_pnd_LEV ) THEN343 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 344 344 ! concentration 345 345 zamsk = 1._wp … … 361 361 ! 362 362 ! --- Lateral boundary conditions --- ! 363 IF ( ln_pnd_LEV.AND. ln_pnd_lids ) THEN363 IF ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 364 364 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 365 365 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 366 ELSEIF( ln_pnd_LEV.AND. .NOT.ln_pnd_lids ) THEN366 ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 367 367 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 368 368 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) … … 1611 1611 ! ! -- check h_ip -- ! 1612 1612 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1613 IF( ln_pnd_LEV.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN1613 IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1614 1614 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1615 1615 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icedyn_rdgrft.F90
r13617 r13902 567 567 oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft 568 568 569 IF ( ln_pnd_LEV ) THEN569 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 570 570 aprdg1 = a_ip_2d(ji,jl1) * afrdg 571 571 aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) … … 604 604 sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1 - sirft(ji) 605 605 oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1 - oirft1 606 IF ( ln_pnd_LEV ) THEN606 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 607 607 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 608 608 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) … … 701 701 v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + & 702 702 & vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 703 IF ( ln_pnd_LEV ) THEN703 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 704 704 v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) & 705 705 & + vprft (ji) * rn_fpndrft * zswitch(ji) ) -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/iceistate.F90
r13284 r13902 381 381 382 382 ! Melt ponds 383 WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 383 WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 384 384 ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp 385 385 END WHERE -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/iceitd.F90
r13617 r13902 305 305 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 306 306 a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 307 IF( ln_pnd_LEV ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin307 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 308 308 h_i_1d(ji) = rn_himin 309 309 ENDIF … … 476 476 zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans 477 477 ! 478 IF ( ln_pnd_LEV ) THEN478 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 479 479 ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction 480 480 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
r13891 r13902 24 24 ! 25 25 USE in_out_manager ! I/O manager 26 USE iom ! I/O manager library 26 27 USE lib_mpp ! MPP library 27 28 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) … … 45 46 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 46 47 zvp_min = 1.e-4_wp ! minimum pond volume (m) 48 49 !-------------------------------------------------------------------------- 47 50 ! 48 51 ! Pond volume per area budget diags … … 58 61 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 59 62 ! 60 ! In level mode 61 62 REAL(wp), DIMENSION(jpi,jpj) ::! pond volume budget diagnostics63 ! In level mode, all terms are incorporated 64 65 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 63 66 diag_dvpn_mlt, & !! meltwater pond volume input [m/day] 64 67 diag_dvpn_drn, & !! pond volume lost by drainage [m/day] 65 diag_dvpn_lid, & !! pond volume lost byrefreezing [m/day]68 diag_dvpn_lid, & !! exchange with lid / refreezing [m/day] 66 69 diag_dvpn_rnf !! meltwater pond lost to runoff [m/day] 67 70 68 REAL(wp), DIMENSION(jpij) ::! pond volume budget diagnostics (1d)71 REAL(wp), ALLOCATABLE, DIMENSION(:) :: & ! pond volume budget diagnostics (1d) 69 72 diag_dvpn_mlt_1d, & !! meltwater pond volume input [m/day] 70 73 diag_dvpn_drn_1d, & !! pond volume lost by drainage [m/day] 71 diag_dvpn_lid_1d, & !! pond volume lost byrefreezing [m/day]74 diag_dvpn_lid_1d, & !! exchange with lid / refreezing [m/day] 72 75 diag_dvpn_rnf_1d !! meltwater pond lost to runoff [m/day] 73 76 ! 77 !-------------------------------------------------------------------------- 74 78 75 79 !! * Substitutions … … 83 87 84 88 SUBROUTINE ice_thd_pnd 89 85 90 !!------------------------------------------------------------------- 86 91 !! *** ROUTINE ice_thd_pnd *** … … 99 104 !! The current diagnostic lacks a contribution from drainage 100 105 !! 101 !! Each time wfx_pnd is updated, wfx_sum / wfx_snw_sum must be updated102 !!103 106 !!------------------------------------------------------------------- 104 ! 105 diag_dvpn_mlt(:,:) = 0._wp ; diag_dvpn_drn(:,:) = 0._wp 106 diag_dvpn_lid(:,:) = 0._wp ; diag_dvpn_rnf(:,:) = 0._wp 107 diag_dvpn_mlt_1d(:,:) = 0._wp ; diag_dvpn_drn_1d(:,:) = 0._wp 108 diag_dvpn_lid_1d(:,:) = 0._wp ; diag_dvpn_rnf_1d(:,:) = 0._wp 107 !! 108 109 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 110 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 111 112 diag_dvpn_mlt(:,:) = 0._wp ; diag_dvpn_drn(:,:) = 0._wp 113 diag_dvpn_lid(:,:) = 0._wp ; diag_dvpn_rnf(:,:) = 0._wp 114 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 115 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 109 116 110 117 SELECT CASE ( nice_pnd ) … … 119 126 ! 120 127 121 IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', dvpn_mlt * 100._wp ) ! input from melting 122 IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', dvpn_lid * 100._wp ) ! exchanges with lid 123 IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', dvpn_drn * 100._wp ) ! vertical drainage 124 IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', dvpn_rnf * 100._wp ) ! runoff + overflow 128 IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 129 IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 130 IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 131 IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 132 133 DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 134 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 125 135 126 136 END SUBROUTINE ice_thd_pnd … … 209 219 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 210 220 !! 211 !! ** Note : mostly stolen from CICE but not only212 !! 213 !! ** References : Holland et al (J. Clim,2012)221 !! ** Note : Mostly stolen from CICE but not only 222 !! 223 !! ** References : Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 214 224 !! 215 225 !!------------------------------------------------------------------- … … 230 240 REAL(wp) :: zfac, zdum ! temporary arrays 231 241 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 232 REAL(wp) :: z1_rdtice ! inverse time step233 242 REAL(wp) :: zvold ! 234 243 !! … … 240 249 z1_aspect = 1._wp / zaspect 241 250 z1_Tp = 1._wp / zTp 242 z1_rdtice = 1._wp / rdt_ice243 251 244 252 !----------------------------------------------------------------------------------------------- … … 291 299 ! Go for ponds 292 300 !----------------------------------------------------------------------------------------------- 293 294 301 295 302 DO ji = 1, npti … … 316 323 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 317 324 318 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * z1_rdtice ! surface melt input diag 319 320 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * z1_rdtice ! runoff diag 321 325 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_rdtice ! surface melt input diag 326 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_rdtice ! runoff diag 322 327 ! 323 328 !--- overflow ---! … … 339 344 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 340 345 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 341 342 346 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 343 344 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * z1_rdtice ! runoff diag - overflow contribution 347 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_rdtice ! runoff diag - overflow contribution 345 348 346 349 !--- Pond growing ---! … … 391 394 ENDIF 392 395 393 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * z1_rdtice ! shrinking counted as lid diagnostic396 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_rdtice ! shrinking counted as lid diagnostic 394 397 395 398 ! … … 419 422 420 423 ! Do the drainage using Darcy's law 421 zdv_flush = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)424 zdv_flush = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 422 425 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 423 426 ! zdv_flush = 0._wp ! MV remove pond drainage for now 424 427 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 425 428 426 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * z1_rdtice ! shrinking counted as lid diagnostic429 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_rdtice ! shrinking counted as lid diagnostic 427 430 428 431 ! MV --- why pond drainage does not give back water into freshwater flux ? … … 486 489 !! *** ROUTINE pnd_TOPO *** 487 490 !! 488 !! ** Purpose : Compute melt pond evolution489 !!490 491 !! ** Purpose : Compute melt pond evolution based on the ice 491 !! topography as inferred from the ice thickness 492 !! distribution. 492 !! topography inferred from the ice thickness distribution 493 493 !! 494 494 !! ** Method : This code is initially based on Flocco and Feltham 495 !! (2007) and Flocco et al. (2010). More to come... 495 !! (2007) and Flocco et al. (2010). 496 !! 497 !! - Calculate available pond water base on surface meltwater 498 !! - Redistribute water as a function of topography, drain water 499 !! - Exchange water with the lid 496 500 !! 497 501 !! ** Tunable parameters : … … 524 528 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 525 529 z1_rhow, & ! inverse water density 526 z1_rdtice, & ! inverse time step527 530 zv_pnd , & ! volume of meltwater contributing to ponds 528 531 zv_mlt ! total amount of meltwater produced … … 550 553 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 551 554 z1_rhow = 1._wp / rhow 552 z1_rdtice = 1._wp / rdt_ice553 555 554 556 ! Set required ice variables (hard-coded here for now) … … 569 571 ! Change 2D to 1D 570 572 !--------------------------------------------------------------- 571 573 ! MV 574 ! a less computing-intensive version would have 2D-1D passage here 572 575 ! use what we have in iceitd.F90 (incremental remapping) 573 576 … … 590 593 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 591 594 592 !--- Available meltwater for melt ponding ---! 595 !--- Available and contributing meltwater for melt ponding ---! 596 zv_mlt = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) & ! available volume of surface melt water per grid area 597 & * z1_rhow * a_i(ji,jj,jl) 598 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 593 599 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 594 595 zv_mlt = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) & ! total volume of surface melt water 596 & * z1_rhow * a_i(ji,jj,jl) 597 zv_pnd = zfr_mlt * zv_mlt ! volume of meltwater contributing to ponds 598 599 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * z1_rdtice ! diagnostics 600 601 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * z1_rdtice 600 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 601 602 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice ! diags 603 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice 602 604 603 605 !--- Create possible new ponds 604 606 ! if pond does not exist, create new pond over full ice area 605 IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 606 h_ip(ji,jj,jl) = 0._wp 607 !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 608 IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 609 a_ip(ji,jj,jl) = a_i(ji,jj,jl) 607 610 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 608 a_ip(ji,jj,jl) = a_i(ji,jj,jl)609 611 ENDIF 610 612 611 !--- Deepen existing ponds before redistribution and drainage 612 ! without pond fraction 613 v_ip(ji,jj,jl) = v_i_p(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 614 613 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 614 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 615 615 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 616 616 617 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl)618 617 !--- Total available pond water volume (pre-existing + newly produced)j 618 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 619 619 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 620 620 … … 631 631 632 632 !-------------------------------------------------------------- 633 ! Freeze and melt lid633 ! Melt pond lid growth and melt 634 634 !-------------------------------------------------------------- 635 DO jj = 1, jpj 636 DO ji = 1, jpi 637 638 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min *at_i(ji,jj) ) THEN 635 636 IF( ln_pnd_lids ) THEN 637 638 DO jj = 1, jpj 639 DO ji = 1, jpi 640 641 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 639 642 640 643 !-------------------------- … … 732 735 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 733 736 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 734 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! in principle, this is useless as h_ip is computed in icevar737 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 735 738 ENDIF 736 739 … … 750 753 END DO ! jj 751 754 END DO ! ji 755 756 ENDIF ! ln_pnd_lids 752 757 753 758 !--------------------------------------------------------------- … … 759 764 760 765 DO jl = 1, jpl 766 761 767 DO jj = 1, jpj 762 768 DO ji = 1, jpi 763 769 764 IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 765 .AND. v_il (ji,jj,jl) > epsi10) THEN 766 v_il(ji,jj,jl) = 0._wp 770 ! ! zap lids on small ponds 771 ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 772 ! .AND. v_il(ji,jj,jl) > epsi10) THEN 773 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 774 ! ENDIF 775 776 ! recalculate equivalent pond variables 777 IF ( a_ip(ji,jj,jl) > epsi10) THEN 778 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 779 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 780 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 767 781 ENDIF 768 769 ! reload tracers 770 IF ( a_ip(ji,jj,jl) > epsi10) THEN 771 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 772 ELSE 773 v_il(ji,jj,jl) = 0._wp 774 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar 775 ENDIF 782 ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 783 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 784 ! ENDIF 776 785 777 IF ( a_ip(ji,jj,jl) > epsi10 ) THEN 778 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 779 !h_ip(ji,jj,jl) = zhpondn(ji,jj,jl) 780 ELSE 781 a_ip_frac(ji,jj,jl) = 0._wp 782 h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 783 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 784 ENDIF 785 786 END DO ! ji 786 END DO ! ji 787 787 END DO ! jj 788 788 789 END DO ! jl 789 790 … … 799 800 !! redistribute and drain water 800 801 !! 801 !! ** 802 !! ** Method 803 !! 804 !-----------| 805 ! | 806 ! |-----------| 807 !___________|___________|______________________________________sea-level 808 ! | | 809 ! | |---^--------| 810 ! | | | | 811 ! | | | |-----------| |------- 812 ! | | | alfan | | | 813 ! | | | | |--------------| 814 ! | | | | | | 815 !---------------------------v------------------------------------------- 816 ! | | ^ | | | 817 ! | | | | |--------------| 818 ! | | | betan | | | 819 ! | | | |-----------| |------- 820 ! | | | | 821 ! | |---v------- | 822 ! | | 823 ! |-----------| 824 ! | 825 !-----------| 826 ! 802 827 !! 803 828 !!------------------------------------------------------------------ 804 829 805 REAL (wp), DIMENSION(jpi,jpj), 806 zvolp, 807 zdvolp 808 809 INTEGER :: &810 ns, &830 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 831 zvolp, & ! total available pond water 832 zdvolp ! remaining meltwater after redistribution 833 834 INTEGER :: & 835 ns, & 811 836 m_index, & 812 837 permflag … … 839 864 INTEGER :: ji, jj, jk, jl ! loop indices 840 865 841 !-----------| 842 ! | 843 ! |-----------| 844 !___________|___________|______________________________________sea-level 845 ! | | 846 ! | |---^--------| 847 ! | | | | 848 ! | | | |-----------| |------- 849 ! | | |alfan(jl)| | | 850 ! | | | | |--------------| 851 ! | | | | | | 852 !---------------------------v------------------------------------------- 853 ! | | ^ | | | 854 ! | | | | |--------------| 855 ! | | |betan(jl)| | | 856 ! | | | |-----------| |------- 857 ! | | | | 858 ! | |---v------- | 859 ! | | 860 ! |-----------| 861 ! | 862 !-----------| 863 864 a_ip(:,:,:) = 0._wp 865 h_ip(:,:,:) = 0._wp 866 867 DO jj = 1, jpj 868 DO ji = 1, jpi 869 870 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 871 872 !------------------------------------------------------------------- 873 ! initialize 874 !------------------------------------------------------------------- 875 876 DO jl = 1, jpl 877 878 a_ip(ji,jj,jl) = 0._wp 879 h_ip(ji,jj,jl) = 0._wp 880 881 !---------------------------------------- 882 ! compute the effective snow fraction 883 !---------------------------------------- 884 885 IF (a_i(ji,jj,jl) < epsi10) THEN 886 hicen(jl) = 0._wp 887 hsnon(jl) = 0._wp 888 reduced_aicen(jl) = 0._wp 889 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 866 a_ip(:,:,:) = 0._wp 867 h_ip(:,:,:) = 0._wp 868 869 DO jj = 1, jpj 870 DO ji = 1, jpi 871 872 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 873 874 !------------------------------------------------------------------- 875 ! initialize 876 !------------------------------------------------------------------- 877 878 DO jl = 1, jpl 879 880 !---------------------------------------- 881 ! compute the effective snow fraction 882 !---------------------------------------- 883 884 IF (a_i(ji,jj,jl) < epsi10) THEN 885 hicen(jl) = 0._wp 886 hsnon(jl) = 0._wp 887 reduced_aicen(jl) = 0._wp 888 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 889 ELSE 890 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 891 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 892 reduced_aicen(jl) = 1._wp ! n=jpl 893 894 !js: initial code in NEMO_DEV 895 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 896 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 897 898 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 899 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 900 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 901 902 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 903 ! MV should check whether this makes sense to have the same effective snow fraction in here 904 ! OLI: it probably doesn't 905 END IF 906 907 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 908 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 909 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 910 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 911 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 912 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 913 ! alfan = 60% of the ice volume) in each category lies above the reference line, 914 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 915 916 ! MV: 917 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 918 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 919 920 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 921 922 alfan(jl) = 0.6 * hicen(jl) 923 betan(jl) = 0.4 * hicen(jl) 924 925 cum_max_vol(jl) = 0._wp 926 cum_max_vol_tmp(jl) = 0._wp 927 928 END DO ! jpl 929 930 cum_max_vol_tmp(0) = 0._wp 931 drain = 0._wp 932 zdvolp(ji,jj) = 0._wp 933 934 !---------------------------------------------------------- 935 ! Drain overflow water, update pond fraction and volume 936 !---------------------------------------------------------- 937 938 !-------------------------------------------------------------------------- 939 ! the maximum amount of water that can be contained up to each ice category 940 !-------------------------------------------------------------------------- 941 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 942 ! Then the excess volume cum_max_vol(jl) drains out of the system 943 ! It should be added to wfx_pnd_out 944 945 DO jl = 1, jpl-1 ! last category can not hold any volume 946 947 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 948 949 ! total volume in level including snow 950 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 951 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 952 953 ! subtract snow solid volumes from lower categories in current level 954 DO ns = 1, jl 955 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 956 - rhos/rhow * & ! free air fraction that can be filled by water 957 asnon(ns) * & ! effective areal fraction of snow in that category 958 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 959 END DO 960 961 ELSE ! assume higher categories unoccupied 962 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 963 END IF 964 !IF (cum_max_vol_tmp(jl) < z0) THEN 965 ! CALL abort_ice('negative melt pond volume') 966 !END IF 967 END DO 968 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 969 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 970 971 !---------------------------------------------------------------- 972 ! is there more meltwater than can be held in the floe? 973 !---------------------------------------------------------------- 974 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 975 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 976 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 977 978 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 979 980 zdvolp(ji,jj) = drain ! this is the drained water 981 IF (zvolp(ji,jj) < epsi10) THEN 982 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 983 zvolp(ji,jj) = 0._wp 984 END IF 985 END IF 986 987 ! height and area corresponding to the remaining volume 988 ! routine leaves zvolp unchanged 989 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 990 991 DO jl = 1, m_index 992 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 993 ! ! volume instead, no ? 994 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 995 a_ip(ji,jj,jl) = reduced_aicen(jl) 996 ! in practise, pond fraction depends on the empirical snow fraction 997 ! so in turn on ice thickness 998 END DO 999 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1000 1001 !------------------------------------------------------------------------ 1002 ! Drainage through brine network (permeability) 1003 !------------------------------------------------------------------------ 1004 !!! drainage due to ice permeability - Darcy's law 1005 1006 ! sea water level 1007 msno = 0._wp 1008 DO jl = 1 , jpl 1009 msno = msno + v_s(ji,jj,jl) * rhos 1010 END DO 1011 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 1012 hsl_rel = floe_weight / rau0 & 1013 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1014 1015 deltah = hpond - hsl_rel 1016 pressure_head = grav * rau0 * max(deltah, 0._wp) 1017 1018 ! drain if ice is permeable 1019 permflag = 0 1020 1021 IF (pressure_head > 0._wp) THEN 1022 DO jl = 1, jpl-1 1023 IF ( hicen(jl) /= 0._wp ) THEN 1024 1025 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1026 1027 perm = 0._wp ! MV ugly dummy patch 1028 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1029 IF (perm > 0._wp) permflag = 1 1030 1031 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 1032 (viscosity*hicen(jl)) 1033 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1034 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1035 1036 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1037 1038 IF (zvolp(ji,jj) < epsi10) THEN 1039 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1040 zvolp(ji,jj) = 0._wp 1041 END IF 1042 END IF 1043 END DO 1044 1045 ! adjust melt pond dimensions 1046 IF (permflag > 0) THEN 1047 ! recompute pond depth 1048 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1049 DO jl = 1, m_index 1050 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1051 a_ip(ji,jj,jl) = reduced_aicen(jl) 1052 END DO 1053 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1054 END IF 1055 END IF ! pressure_head 1056 1057 !------------------------------- 1058 ! remove water from the snow 1059 !------------------------------- 1060 !------------------------------------------------------------------------ 1061 ! total melt pond volume in category does not include snow volume 1062 ! snow in melt ponds is not melted 1063 !------------------------------------------------------------------------ 1064 1065 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1066 ! how much, so I did not diagnose it 1067 ! so if there is a problem here, nobody is going to see it... 1068 1069 1070 ! Calculate pond volume for lower categories 1071 DO jl = 1,m_index-1 1072 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1073 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1074 END DO 1075 1076 ! Calculate pond volume for highest category = remaining pond volume 1077 1078 ! The following is completely unclear to Martin at least 1079 ! Could we redefine properly and recode in a more readable way ? 1080 1081 ! m_index = last category with melt pond 1082 1083 IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 1084 1085 IF (m_index > 1) THEN 1086 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1087 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 890 1088 ELSE 891 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 892 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 893 reduced_aicen(jl) = 1._wp ! n=jpl 894 895 !js: initial code in NEMO_DEV 896 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 897 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 898 899 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 900 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 901 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 902 903 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 904 ! MV should check whether this makes sense to have the same effective snow fraction in here 905 ! OLI: it probably doesn't 1089 v_ip(ji,jj,m_index) = 0._wp 1090 h_ip(ji,jj,m_index) = 0._wp 1091 a_ip(ji,jj,m_index) = 0._wp 1092 ! If remaining pond volume is negative reduce pond volume of 1093 ! lower category 1094 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1095 v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 906 1096 END IF 907 908 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 909 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 910 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 911 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 912 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 913 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 914 ! alfan = 60% of the ice volume) in each category lies above the reference line, 915 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 916 917 ! MV: 918 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 919 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 920 921 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 922 923 alfan(jl) = 0.6 * hicen(jl) 924 betan(jl) = 0.4 * hicen(jl) 925 926 cum_max_vol(jl) = 0._wp 927 cum_max_vol_tmp(jl) = 0._wp 928 929 END DO ! jpl 930 931 cum_max_vol_tmp(0) = 0._wp 932 drain = 0._wp 933 zdvolp(ji,jj) = 0._wp 934 935 !---------------------------------------------------------- 936 ! Drain overflow water, update pond fraction and volume 937 !---------------------------------------------------------- 938 939 !-------------------------------------------------------------------------- 940 ! the maximum amount of water that can be contained up to each ice category 941 !-------------------------------------------------------------------------- 942 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 943 ! Then the excess volume cum_max_vol(jl) drains out of the system 944 ! It should be added to wfx_pnd_out 945 946 DO jl = 1, jpl-1 ! last category can not hold any volume 947 948 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 949 950 ! total volume in level including snow 951 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 952 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 953 954 ! subtract snow solid volumes from lower categories in current level 955 DO ns = 1, jl 956 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 957 - rhos/rhow * & ! free air fraction that can be filled by water 958 asnon(ns) * & ! effective areal fraction of snow in that category 959 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 960 END DO 961 962 ELSE ! assume higher categories unoccupied 963 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 964 END IF 965 !IF (cum_max_vol_tmp(jl) < z0) THEN 966 ! CALL abort_ice('negative melt pond volume') 967 !END IF 968 END DO 969 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 970 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 971 972 !---------------------------------------------------------------- 973 ! is there more meltwater than can be held in the floe? 974 !---------------------------------------------------------------- 975 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 976 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 977 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 978 979 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 980 981 zdvolp(ji,jj) = drain ! this is the drained water 982 IF (zvolp(ji,jj) < epsi10) THEN 983 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 984 zvolp(ji,jj) = 0._wp 985 END IF 986 END IF 987 988 ! height and area corresponding to the remaining volume 989 ! routine leaves zvolp unchanged 990 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 991 992 DO jl = 1, m_index 993 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 994 ! ! volume instead, no ? 995 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 996 a_ip(ji,jj,jl) = reduced_aicen(jl) 997 ! in practise, pond fraction depends on the empirical snow fraction 998 ! so in turn on ice thickness 999 END DO 1000 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1001 1002 !------------------------------------------------------------------------ 1003 ! Drainage through brine network (permeability) 1004 !------------------------------------------------------------------------ 1005 !!! drainage due to ice permeability - Darcy's law 1006 1007 ! sea water level 1008 msno = 0._wp 1009 DO jl = 1 , jpl 1010 msno = msno + v_s(ji,jj,jl) * rhos 1011 END DO 1012 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 1013 hsl_rel = floe_weight / rau0 & 1014 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1015 1016 deltah = hpond - hsl_rel 1017 pressure_head = grav * rau0 * max(deltah, 0._wp) 1018 1019 ! drain if ice is permeable 1020 permflag = 0 1021 1022 IF (pressure_head > 0._wp) THEN 1023 DO jl = 1, jpl-1 1024 IF ( hicen(jl) /= 0._wp ) THEN 1025 1026 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1027 1028 perm = 0._wp ! MV ugly dummy patch 1029 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1030 IF (perm > 0._wp) permflag = 1 1031 1032 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 1033 (viscosity*hicen(jl)) 1034 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1035 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1036 1037 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1038 1039 IF (zvolp(ji,jj) < epsi10) THEN 1040 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1041 zvolp(ji,jj) = 0._wp 1042 END IF 1043 END IF 1044 END DO 1045 1046 ! adjust melt pond dimensions 1047 IF (permflag > 0) THEN 1048 ! recompute pond depth 1049 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1050 DO jl = 1, m_index 1051 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1052 a_ip(ji,jj,jl) = reduced_aicen(jl) 1053 END DO 1054 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1055 END IF 1056 END IF ! pressure_head 1057 1058 !------------------------------- 1059 ! remove water from the snow 1060 !------------------------------- 1061 !------------------------------------------------------------------------ 1062 ! total melt pond volume in category does not include snow volume 1063 ! snow in melt ponds is not melted 1064 !------------------------------------------------------------------------ 1065 1066 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1067 ! how much, so I did not diagnose it 1068 ! so if there is a problem here, nobody is going to see it... 1069 1070 1071 ! Calculate pond volume for lower categories 1072 DO jl = 1,m_index-1 1073 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1074 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1075 END DO 1076 1077 ! Calculate pond volume for highest category = remaining pond volume 1078 1079 ! The following is completely unclear to Martin at least 1080 ! Could we redefine properly and recode in a more readable way ? 1081 1082 ! m_index = last category with melt pond 1083 1084 IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 1085 1086 IF (m_index > 1) THEN 1087 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1088 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1089 ELSE 1090 v_ip(ji,jj,m_index) = 0._wp 1091 h_ip(ji,jj,m_index) = 0._wp 1092 a_ip(ji,jj,m_index) = 0._wp 1093 ! If remaining pond volume is negative reduce pond volume of 1094 ! lower category 1095 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1096 v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 1097 END IF 1098 END IF 1099 1100 DO jl = 1,m_index 1101 IF (a_ip(ji,jj,jl) > epsi10) THEN 1102 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1103 ELSE 1104 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1105 h_ip(ji,jj,jl) = 0._wp 1106 v_ip(ji,jj,jl) = 0._wp 1107 a_ip(ji,jj,jl) = 0._wp 1108 END IF 1109 END DO 1110 DO jl = m_index+1, jpl 1111 h_ip(ji,jj,jl) = 0._wp 1112 a_ip(ji,jj,jl) = 0._wp 1113 v_ip(ji,jj,jl) = 0._wp 1114 END DO 1115 1116 ENDIF 1117 END DO ! ji 1118 END DO ! jj 1097 END IF 1098 1099 DO jl = 1,m_index 1100 IF (a_ip(ji,jj,jl) > epsi10) THEN 1101 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1102 ELSE 1103 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1104 h_ip(ji,jj,jl) = 0._wp 1105 v_ip(ji,jj,jl) = 0._wp 1106 a_ip(ji,jj,jl) = 0._wp 1107 END IF 1108 END DO 1109 DO jl = m_index+1, jpl 1110 h_ip(ji,jj,jl) = 0._wp 1111 a_ip(ji,jj,jl) = 0._wp 1112 v_ip(ji,jj,jl) = 0._wp 1113 END DO 1114 1115 ENDIF 1116 END DO ! ji 1117 END DO ! jj 1119 1118 1120 1119 END SUBROUTINE ice_thd_pnd_area … … 1355 1354 INTEGER :: ios, ioptio ! Local integer 1356 1355 !! 1357 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &1358 & ln_pnd_CST , rn_apnd, rn_hpnd, &1359 & ln_pnd_TOPO , &1356 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 1357 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1358 & ln_pnd_TOPO , & 1360 1359 & ln_pnd_lids, ln_pnd_alb 1361 1360 !!------------------------------------------------------------------- … … 1379 1378 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 1380 1379 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 1380 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush 1381 1381 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 1382 1382 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd … … 1393 1393 IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF 1394 1394 IF( ioptio /= 1 ) & 1395 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )1395 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' ) 1396 1396 ! 1397 1397 SELECT CASE( nice_pnd ) -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icevar.F90
r13284 r13902 565 565 END DO 566 566 END DO 567 568 !----------------------------------------------------------------- 569 ! zap small ponds 570 !----------------------------------------------------------------- 571 DO jj = 1, jpj 572 DO ji = 1, jpi 573 IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN 574 a_ip(ji,jj,jl) = 0._wp 575 a_ip_frac(ji,jj,jl) = 0._wp 576 v_ip(ji,jj,jl) = 0._wp 577 h_ip(ji,jj,jl) = 0._wp 578 v_il(ji,jj,jl) = 0._wp 579 h_il(ji,jj,jl) = 0._wp 580 ENDIF 581 END DO 582 END DO 567 583 ! 568 584 END DO … … 696 712 WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 697 713 WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 698 IF( ln_pnd_LEV ) THEN714 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 699 715 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 700 716 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icewri.F90
r13284 r13902 164 164 IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 165 165 IF( iom_use('iceapnd_cat' ) ) CALL iom_put( 'iceapnd_cat' , a_ip * zmsk00l ) ! melt pond frac for categories 166 IF( iom_use('icevpnd_cat' ) ) CALL iom_put( 'icevpnd_cat' , v_ip * zmsk00l ) ! melt pond volume for categories 166 167 IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 167 168 IF( iom_use('icehlid_cat' ) ) CALL iom_put( 'icehlid_cat' , h_il * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 168 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac for categories169 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac per ice area for categories 169 170 IF( iom_use('iceaepnd_cat') ) CALL iom_put( 'iceaepnd_cat', a_ip_eff * zmsk00l ) ! melt pond effective frac for categories 170 171 IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
Note: See TracChangeset
for help on using the changeset viewer.
