Changeset 13957 for NEMO/branches/2020/SI3_martin_ponds/src
- Timestamp:
- 2020-12-01T23:04:18+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/SI3_martin_ponds/src/ICE
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90
r13908 r13957 311 311 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: cnd_ice !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 312 312 313 ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array)314 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_i_sum_2d !: surface melt (2d arrays for ponds) [m]315 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_s_mlt_2d !: snow surf melt (2d arrays for ponds) [m]316 317 313 !!---------------------------------------------------------------------- 318 314 !! * Ice global state variables … … 368 364 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_il !: total melt pond lid volume per gridcell area [m] 369 365 366 ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array) 367 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_i_sum_2d !: surface melt (2d arrays for ponds) [m] 368 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_s_mlt_2d !: snow surf melt (2d arrays for ponds) [m] 369 370 370 !!---------------------------------------------------------------------- 371 371 !! * Global variables at before time step 372 372 !!---------------------------------------------------------------------- 373 373 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 374 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip_b, v_il_b !: ponds and lids volumes 374 375 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b !: 375 376 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content … … 398 399 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 399 400 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 401 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vpnd !: pond volume variation [m/s] 400 402 ! 401 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) … … 464 466 ii = ii + 1 465 467 ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , & 466 & dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , &467 468 & h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , & 468 469 & v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & … … 485 486 ii = ii + 1 486 487 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), & 487 & v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 488 & v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , & 489 & dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) ) 488 490 489 491 ii = ii + 1 … … 493 495 ii = ii + 1 494 496 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), & 497 & v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , & 495 498 & a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 496 499 & STAT=ierr(ii) ) … … 507 510 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 508 511 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 509 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), &512 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 510 513 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 511 514 -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icectl.F90
r13601 r13957 85 85 !! 86 86 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, & 87 & zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 87 & zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 88 & zdiag_eimin, zdiag_esmin, zdiag_simin 88 89 REAL(wp) :: zvtrp, zetrp 89 90 REAL(wp) :: zarea … … 92 93 IF( icount == 0 ) THEN 93 94 94 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos , dim=3 ) * e1e2t )95 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 95 96 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 96 97 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) … … 112 113 113 114 ! -- mass diag -- ! 114 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_Dt_ice & 115 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) & 116 & - pdiag_v ) * r1_Dt_ice & 115 117 & + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + & 116 118 & wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & … … 132 134 133 135 ! -- min/max diag -- ! 134 zdiag_amax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 135 zdiag_vmin = glob_min( 'icectl', v_i ) 136 zdiag_amin = glob_min( 'icectl', a_i ) 137 zdiag_smin = glob_min( 'icectl', sv_i ) 136 zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 137 zdiag_vimin = glob_min( 'icectl', v_i ) 138 zdiag_vsmin = glob_min( 'icectl', v_s ) 139 zdiag_vpmin = glob_min( 'icectl', v_ip ) 140 zdiag_vlmin = glob_min( 'icectl', v_il ) 141 zdiag_aimin = glob_min( 'icectl', a_i ) 142 zdiag_simin = glob_min( 'icectl', sv_i ) 138 143 zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 139 144 zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) … … 155 160 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 156 161 ! check negative values 157 IF( zdiag_vmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vmin 158 IF( zdiag_amin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_amin 159 IF( zdiag_smin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_smin 160 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 161 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 IF( zdiag_vimin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vimin 163 IF( zdiag_vsmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_s < 0 = ',zdiag_vsmin 164 IF( zdiag_vpmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_ip < 0 = ',zdiag_vpmin 165 IF( zdiag_vlmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_il < 0 = ',zdiag_vlmin 166 IF( zdiag_aimin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_aimin 167 IF( zdiag_simin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_simin 168 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 169 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 170 ! check maximum ice concentration 163 IF( zdiag_a max >MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_a max171 IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 172 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_aimax 165 173 ! check if advection scheme is conservative 166 174 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * r dt_ice175 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 168 176 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * r dt_ice177 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rDt_ice 170 178 ENDIF 171 179 ! … … 193 201 ! water flux 194 202 ! -- mass diag -- ! 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub&196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t )203 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd & 204 & + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 197 205 198 206 ! -- salt diag -- ! … … 200 208 201 209 ! -- heat diag -- ! 202 zdiag_heat 210 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 211 ! equivalent to this: 204 212 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & … … 245 253 IF( icount == 0 ) THEN 246 254 247 pdiag_v = SUM( v_i * rhoi + v_s * rhos , dim=3 )248 pdiag_s = SUM( sv_i * rhoi 255 pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 256 pdiag_s = SUM( sv_i * rhoi , dim=3 ) 249 257 pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 250 258 … … 261 269 262 270 ! -- mass diag -- ! 263 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos , dim=3 ) - pdiag_v ) * r1_Dt_ice&271 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice & 264 272 & + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 265 273 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & … … 352 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 353 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 362 ! mean state 363 CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) ! 364 CALL iom_rstput( 0, 0, inum, 'icevol' , SUM(v_i ,dim=3) , ktype = jp_r8 ) ! 365 CALL iom_rstput( 0, 0, inum, 'snwvol' , SUM(v_s ,dim=3) , ktype = jp_r8 ) ! 366 CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) ! 367 CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) ! 354 368 355 369 CALL iom_close( inum ) … … 776 790 ! -- mass diag -- ! 777 791 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * r dt_ice792 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rDt_ice 779 793 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 794 781 795 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * r dt_ice * 1.e-3796 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 783 797 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 798 -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_pra.F90
r13908 r13957 156 156 157 157 ! diagnostics 158 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 158 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 159 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 159 160 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 160 161 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 320 321 ! 321 322 ! --- diagnostics --- ! 322 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 323 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 324 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 323 325 & - zdiag_adv_mass(:,:) ) * z1_dt 324 326 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90
r13911 r13957 182 182 183 183 ! diagnostics 184 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 184 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 185 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 185 186 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 186 187 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 379 380 ! 380 381 ! --- diagnostics --- ! 381 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 382 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 383 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 382 384 & - zdiag_adv_mass(:,:) ) * z1_dt 383 385 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceistate.F90
r13472 r13957 426 426 ! 4) Snow-ice mass (case ice is fully embedded) 427 427 !---------------------------------------------- 428 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s (:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass428 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3 ) ! snow+ice mass 429 429 snwice_mass_b(:,:) = snwice_mass(:,:) 430 430 ! -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceitd.F90
r13908 r13957 29 29 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 30 30 USE prtctl ! Print control 31 USE timing ! Timing 31 32 32 33 IMPLICIT NONE … … 87 88 REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories 88 89 !!------------------------------------------------------------------ 90 IF( ln_timing ) CALL timing_start('iceitd_rem') 89 91 90 92 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' … … 328 330 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 329 331 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 332 IF( ln_timing ) CALL timing_stop ('iceitd_rem') 330 333 ! 331 334 END SUBROUTINE ice_itd_rem … … 491 494 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 492 495 ! 493 ztrans = v_ip_2d(ji,jl1) * zwork a(ji) ! Pond volume (also proportional to da/a)496 ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume 494 497 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 495 498 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 496 499 ! 497 500 IF ( ln_pnd_lids ) THEN ! Pond lid volume 498 ztrans = v_il_2d(ji,jl1) * zwork a(ji)501 ztrans = v_il_2d(ji,jl1) * zworkv(ji) 499 502 v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 500 503 v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans … … 606 609 REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred 607 610 !!------------------------------------------------------------------ 611 IF( ln_timing ) CALL timing_start('iceitd_reb') 608 612 ! 609 613 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' … … 635 639 jdonor(ji,jl) = jl 636 640 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 641 ! these are from CICE => transfer everything 642 !!zdaice(ji,jl) = a_i_1d(ji) 643 !!zdvice(ji,jl) = v_i_1d(ji) 644 ! these are from LLN => transfer only half of the category 645 zdaice(ji,jl) = 0.5_wp * a_i_1d(ji) 646 zdvice(ji,jl) = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 646 647 END DO 647 648 ! … … 686 687 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 687 688 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 689 IF( ln_timing ) CALL timing_stop ('iceitd_reb') 688 690 ! 689 691 END SUBROUTINE ice_itd_reb -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icesbc.F90
r13472 r13957 62 62 !!------------------------------------------------------------------- 63 63 ! 64 IF( ln_timing ) CALL timing_start('ice _sbc')64 IF( ln_timing ) CALL timing_start('icesbc') 65 65 ! 66 66 IF( kt == nit000 .AND. lwp ) THEN … … 89 89 ENDIF 90 90 ! 91 IF( ln_timing ) CALL timing_stop('ice _sbc')91 IF( ln_timing ) CALL timing_stop('icesbc') 92 92 ! 93 93 END SUBROUTINE ice_sbc_tau … … 122 122 !!-------------------------------------------------------------------- 123 123 ! 124 IF( ln_timing ) CALL timing_start('ice _sbc_flx')124 IF( ln_timing ) CALL timing_start('icesbc') 125 125 126 126 IF( kt == nit000 .AND. lwp ) THEN … … 176 176 ENDIF 177 177 ! 178 IF( ln_timing ) CALL timing_stop('ice _sbc_flx')178 IF( ln_timing ) CALL timing_stop('icesbc') 179 179 ! 180 180 END SUBROUTINE ice_sbc_flx -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icestp.F90
r13908 r13957 121 121 !!---------------------------------------------------------------------- 122 122 ! 123 IF( ln_timing ) CALL timing_start('ice _stp')123 IF( ln_timing ) CALL timing_start('icestp') 124 124 ! 125 125 ! !-----------------------! … … 215 215 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 216 216 ! 217 IF( ln_timing ) CALL timing_stop('ice _stp')217 IF( ln_timing ) CALL timing_stop('icestp') 218 218 ! 219 219 END SUBROUTINE ice_stp … … 370 370 v_i_b (:,:,:) = v_i (:,:,:) ! ice volume 371 371 v_s_b (:,:,:) = v_s (:,:,:) ! snow volume 372 v_ip_b(:,:,:) = v_ip(:,:,:) ! pond volume 373 v_il_b(:,:,:) = v_il(:,:,:) ! pond lid volume 372 374 sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content 373 375 e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy … … 429 431 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 430 432 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 433 diag_aice(ji,jj) = 0._wp ; diag_vpnd(ji,jj) = 0._wp 431 434 432 435 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) … … 485 488 diag_vsnw(:,:) = diag_vsnw(:,:) & 486 489 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos 490 diag_vpnd(:,:) = diag_vpnd(:,:) & 491 & + SUM( v_ip + v_il - v_ip_b - v_il_b , dim=3 ) * r1_Dt_ice * rhow 487 492 ! 488 493 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90
r13908 r13957 259 259 IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 260 260 ! 261 IF ( ln_pnd ) CALL ice_thd_pnd ! --- Melt ponds 261 IF ( ln_pnd .AND. ln_icedH ) & 262 & CALL ice_thd_pnd ! --- Melt ponds 262 263 ! 263 264 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! … … 267 268 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 268 269 ! 269 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * r dt_ice ! ice natural aging incrementation270 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation 270 271 ! 271 272 ! convergence tests … … 378 379 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 379 380 END DO 380 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )381 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )382 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )383 381 ! 384 382 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 410 408 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 411 409 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 412 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )413 410 ! 414 411 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 465 462 v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 466 463 sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 467 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)468 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)469 464 oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 470 465 … … 484 479 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 485 480 END DO 486 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )487 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )488 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )489 481 ! 490 482 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 502 494 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 503 495 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 504 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )505 496 ! 506 497 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 541 532 CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 542 533 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 543 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )544 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )545 534 CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 546 535 ! check convergence of heat diffusion scheme -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90
r13931 r13957 36 36 INTEGER :: nice_pnd ! choice of the type of pond scheme 37 37 ! ! associated indices: 38 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme39 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme38 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 39 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 41 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 42 42 43 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 44 zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m) 45 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 46 zvp_min = 1.e-4_wp ! minimum pond volume (m) 47 48 !-------------------------------------------------------------------------- 49 ! 50 ! Pond volume per area budget diags 51 ! 52 ! The idea of diags is the volume of ponds per grid cell area is 43 !-------------------------------------------------------------------------- 44 ! Diagnostics for pond volume per area 53 45 ! 54 46 ! dV/dt = mlt + drn + lid + rnf … … 59 51 ! 60 52 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 53 ! In level mode, all terms are incorporated 61 54 ! 62 ! In level mode, all terms are incorporated 63 64 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 65 diag_dvpn_mlt, & !! meltwater pond volume input [m/day] 66 diag_dvpn_drn, & !! pond volume lost by drainage [m/day] 67 diag_dvpn_lid, & !! exchange with lid / refreezing [m/day] 68 diag_dvpn_rnf !! meltwater pond lost to runoff [m/day] 69 70 REAL(wp), ALLOCATABLE, DIMENSION(:) :: & ! pond volume budget diagnostics (1d) 71 diag_dvpn_mlt_1d, & !! meltwater pond volume input [m/day] 72 diag_dvpn_drn_1d, & !! pond volume lost by drainage [m/day] 73 diag_dvpn_lid_1d, & !! exchange with lid / refreezing [m/day] 74 diag_dvpn_rnf_1d !! meltwater pond lost to runoff [m/day] 55 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_mlt ! meltwater pond volume input [kg/m2/s] 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_drn ! pond volume lost by drainage [-] 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_lid ! exchange with lid / refreezing [-] 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] 59 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_mlt_1d ! meltwater pond volume input [-] 60 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_drn_1d ! pond volume lost by drainage [-] 61 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_lid_1d ! exchange with lid / refreezing [-] 62 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_rnf_1d ! meltwater pond lost to runoff [-] 75 63 76 64 !! * Substitutions … … 90 78 !! ** Purpose : change melt pond fraction and thickness 91 79 !! 92 !! Note: Melt ponds affect only radiative transfer for now 93 !! 94 !! No heat, no salt. 95 !! The melt water they carry is collected but 96 !! not removed from fw budget or released to the ocean 97 !! 98 !! A wfx_pnd has been coded for diagnostic purposes 99 !! It is not fully consistent yet. 100 !! 101 !! The current diagnostic lacks a contribution from drainage 102 !! 80 !! ** Note : Melt ponds affect only radiative transfer for now 81 !! No heat, no salt. 82 !! The current diagnostics lacks a contribution from drainage 103 83 !!------------------------------------------------------------------- 104 !! 84 INTEGER :: ji, jj, jl ! loop indices 85 !!------------------------------------------------------------------- 105 86 106 87 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 107 88 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 108 109 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn(:,:)= 0._wp110 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf(:,:)= 0._wp89 ! 90 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn (:,:) = 0._wp 91 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf (:,:) = 0._wp 111 92 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 112 93 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 113 94 114 SELECT CASE ( nice_pnd ) 95 !------------------------------------- 96 ! Remove ponds where ice has vanished 97 !------------------------------------- 98 at_i(:,:) = SUM( a_i, dim=3 ) 115 99 ! 116 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 117 ! 118 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 119 ! 120 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 121 ! 122 END SELECT 100 DO jl = 1, jpl 101 DO_2D( 1, 1, 1, 1 ) 102 IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN 103 wfx_pnd (ji,jj) = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 104 a_ip (ji,jj,jl) = 0._wp 105 v_ip (ji,jj,jl) = 0._wp 106 v_il (ji,jj,jl) = 0._wp 107 h_ip (ji,jj,jl) = 0._wp 108 h_il (ji,jj,jl) = 0._wp 109 a_ip_frac(ji,jj,jl) = 0._wp 110 ENDIF 111 END_2D 112 END DO 113 114 !------------------------------ 115 ! Identify grid cells with ice 116 !------------------------------ 117 npti = 0 ; nptidx(:) = 0 118 DO_2D( 1, 1, 1, 1 ) 119 IF( at_i(ji,jj) >= epsi10 ) THEN 120 npti = npti + 1 121 nptidx( npti ) = (jj - 1) * jpi + ji 122 ENDIF 123 END_2D 124 125 !------------------------------------ 126 ! Select melt pond scheme to be used 127 !------------------------------------ 128 IF( npti > 0 ) THEN 129 SELECT CASE ( nice_pnd ) 130 ! 131 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 132 ! 133 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 134 ! 135 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 136 ! 137 END SELECT 138 ENDIF 139 140 !------------------------------------ 141 ! Diagnostics 142 !------------------------------------ 143 CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting 144 CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid 145 CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage 146 CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow 123 147 ! 124 125 IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 126 IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 127 IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 128 IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 129 130 DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 148 DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf ) 131 149 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 132 150 133 151 END SUBROUTINE ice_thd_pnd 134 135 152 136 153 … … 151 168 !! ** References : Bush, G.W., and Trump, D.J. (2017) 152 169 !!------------------------------------------------------------------- 153 INTEGER :: ji ! loop indices 170 INTEGER :: ji, jl ! loop indices 171 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 154 172 !!------------------------------------------------------------------- 155 DO ji = 1, npti 156 ! 157 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 158 h_ip_1d(ji) = rn_hpnd 159 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 160 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 161 ELSE 162 h_ip_1d(ji) = 0._wp 163 a_ip_1d(ji) = 0._wp 164 h_il_1d(ji) = 0._wp 165 ENDIF 166 ! 173 DO jl = 1, jpl 174 175 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 176 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 177 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 179 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 180 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 181 182 DO ji = 1, npti 183 ! 184 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 185 ! 186 IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 187 h_ip_1d(ji) = rn_hpnd 188 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 189 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 190 ELSE 191 h_ip_1d(ji) = 0._wp 192 a_ip_1d(ji) = 0._wp 193 h_il_1d(ji) = 0._wp 194 ENDIF 195 ! 196 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 197 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 198 ! 199 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 200 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 201 ! 202 END DO 203 204 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 205 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 206 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 207 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 208 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 209 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 210 167 211 END DO 168 212 ! … … 203 247 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 204 248 !! 205 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi 206 !! perm = permability of sea-ice 249 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi * flush --- from Flocco et al 2007 --- 250 !! perm = permability of sea-ice + correction from Hunke et al 2012 (flush) 207 251 !! visc = water viscosity 208 252 !! Hp = height of top of the pond above sea-level 209 253 !! Hi = ice thickness thru which there is flushing 254 !! flush= correction otherwise flushing is excessive 210 255 !! 211 256 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness … … 218 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 219 264 !! 220 !! ** References : Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 221 !! 222 !!------------------------------------------------------------------- 223 265 !! ** References : Flocco and Feltham (JGR, 2007) 266 !! Flocco et al (JGR, 2010) 267 !! Holland et al (J. Clim, 2012) 268 !! Hunke et al (OM 2012) 269 !!------------------------------------------------------------------- 224 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 225 271 !! … … 228 274 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 229 275 !! 230 REAL(wp) :: zfr_mlt, zdv_mlt 276 REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding 231 277 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 278 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 232 279 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 233 280 REAL(wp) :: zv_ip_max ! max pond volume allowed 234 281 REAL(wp) :: zdT ! zTp-t_su 235 REAL(wp) :: zsbr 282 REAL(wp) :: zsbr, ztmelts ! Brine salinity 236 283 REAL(wp) :: zperm ! permeability of sea ice 237 284 REAL(wp) :: zfac, zdum ! temporary arrays 238 285 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 239 REAL(wp) :: zvold ! 240 !! 241 INTEGER :: ji, jj, jk, jl ! loop indices 242 286 !! 287 INTEGER :: ji, jk, jl ! loop indices 243 288 !!------------------------------------------------------------------- 244 245 289 z1_rhow = 1._wp / rhow 246 290 z1_aspect = 1._wp / zaspect 247 291 z1_Tp = 1._wp / zTp 248 292 249 !----------------------------------------------------------------------------------------------- 250 ! Identify grid cells with ice 251 !----------------------------------------------------------------------------------------------- 252 at_i(:,:) = SUM( a_i, dim=3 ) 293 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti), at_i ) 294 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 295 296 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 297 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 298 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 299 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 300 301 DO jl = 1, jpl 302 303 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 304 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 305 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) ) 306 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) ) 307 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) ) 308 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) ) 309 310 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) ) 311 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) ) 312 313 DO jk = 1, nlay_i 314 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 315 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 316 END DO 317 318 !----------------------- 319 ! Melt pond calculations 320 !----------------------- 321 DO ji = 1, npti 322 ! 323 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 324 ! !----------------------------------------------------! 325 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 326 ! !----------------------------------------------------! 327 !--- Remove ponds on thin ice or tiny ice fractions 328 a_ip_1d(ji) = 0._wp 329 h_ip_1d(ji) = 0._wp 330 h_il_1d(ji) = 0._wp 331 ! !--------------------------------! 332 ELSE ! Case ice thickness >= rn_himin ! 333 ! !--------------------------------! 334 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 335 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 336 ! 337 !------------------! 338 ! case ice melting ! 339 !------------------! 340 ! 341 !--- available meltwater for melt ponding (zdv_avail) ---! 342 zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 343 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed 344 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 345 ! 346 !--- overflow ---! 347 ! 348 ! area driven overflow 349 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 350 ! a_ip_max = zfr_mlt * a_i 351 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 352 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 353 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 354 355 ! depth driven overflow 356 ! If pond depth exceeds half the ice thickness then reduce the pond volume 357 ! h_ip_max = 0.5 * h_i 358 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 359 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 360 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 361 362 !--- Pond growing ---! 363 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 364 ! 365 !--- Lid melting ---! 366 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 367 ! 368 !-------------------! 369 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 370 !-------------------! 371 ! 372 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 373 ! 374 !--- Pond contraction (due to refreezing) ---! 375 IF( ln_pnd_lids ) THEN 376 ! 377 !--- Lid growing and subsequent pond shrinking ---! 378 zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 379 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 380 381 ! Lid growing 382 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 383 384 ! Pond shrinking 385 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 386 387 ELSE 388 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 389 ! Pond shrinking 390 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 391 ENDIF 392 ! 393 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 394 ! v_ip = h_ip * a_ip 395 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 396 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 397 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 398 ! 399 400 !------------------------------------------------! 401 ! Pond drainage through brine network (flushing) ! 402 !------------------------------------------------! 403 ! height of top of the pond above sea-level 404 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 405 406 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 407 DO jk = 1, nlay_i 408 ! MV Assur is inconsistent with SI3 409 !!zsbr = - 1.2_wp & 410 !! & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 411 !! & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 412 !! & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 413 !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 414 ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 415 ztmelts = -rTmlt * sz_i_1d(ji,jk) 416 ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 417 END DO 418 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 419 420 ! Do the drainage using Darcy's law 421 zdv_flush = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 423 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 424 425 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 426 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 427 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 428 429 !--- Corrections and lid thickness ---! 430 IF( ln_pnd_lids ) THEN 431 !--- retrieve lid thickness from volume ---! 432 IF( a_ip_1d(ji) > 0.01_wp ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 433 ELSE ; h_il_1d(ji) = 0._wp 434 ENDIF 435 !--- remove ponds if lids are much larger than ponds ---! 436 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 437 a_ip_1d(ji) = 0._wp 438 h_ip_1d(ji) = 0._wp 439 h_il_1d(ji) = 0._wp 440 ENDIF 441 ENDIF 442 443 ! diagnostics: dvpnd = mlt+rnf+lid+drn 444 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow * zdv_avail * r1_Dt_ice ! > 0, surface melt input 445 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice ! < 0, runoff 446 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow * zdv_frz * r1_Dt_ice ! < 0, shrinking 447 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow * zdv_flush * r1_Dt_ice ! < 0, drainage 448 ! 449 ENDIF 450 ! 451 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 452 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 453 ! 454 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 455 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 456 ! 457 END DO 458 459 !-------------------------------------------------------------------- 460 ! Retrieve 2D arrays 461 !-------------------------------------------------------------------- 462 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 466 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 467 ! 468 END DO 253 469 ! 254 npti = 0 ; nptidx(:) = 0 255 DO_2D( 1, 1, 1, 1 ) 256 IF ( at_i(ji,jj) > epsi10 ) THEN 257 npti = npti + 1 258 nptidx( npti ) = (jj - 1) * jpi + ji 259 ENDIF 260 END_2D 261 262 !----------------------------------------------------------------------------------------------- 263 ! Prepare 1D arrays 264 !----------------------------------------------------------------------------------------------- 265 266 IF( npti > 0 ) THEN 267 268 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 269 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti) , wfx_sum ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti) , wfx_pnd ) 271 272 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti) , at_i ) 273 274 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 275 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 276 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 277 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 278 279 DO jl = 1, jpl 280 281 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 282 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 283 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 284 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 285 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 286 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 287 288 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) ) 289 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) ) 290 291 DO jk = 1, nlay_i 292 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 293 END DO 294 295 !----------------------------------------------------------------------------------------------- 296 ! Go for ponds 297 !----------------------------------------------------------------------------------------------- 298 299 DO ji = 1, npti 300 ! !----------------------------------------------------! 301 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 302 ! !----------------------------------------------------! 303 !--- Remove ponds on thin ice or tiny ice fractions 304 a_ip_1d(ji) = 0._wp 305 h_ip_1d(ji) = 0._wp 306 h_il_1d(ji) = 0._wp 307 ! !--------------------------------! 308 ELSE ! Case ice thickness >= rn_himin ! 309 ! !--------------------------------! 310 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 311 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 312 ! 313 !------------------! 314 ! case ice melting ! 315 !------------------! 316 ! 317 !--- available meltwater for melt ponding ---! 318 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 319 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed 320 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 321 322 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_Dt_ice ! surface melt input diag 323 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_Dt_ice ! runoff diag 324 ! 325 !--- overflow ---! 326 ! 327 ! 1) area driven overflow 328 ! 329 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 330 ! a_ip_max = zfr_mlt * a_i 331 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 332 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 333 zvold = zdv_mlt 334 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 335 336 ! 337 ! 2) depth driven overflow 338 ! 339 ! If pond depth exceeds half the ice thickness then reduce the pond volume 340 ! h_ip_max = 0.5 * h_i 341 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 342 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 343 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 344 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_Dt_ice ! runoff diag - overflow contribution 345 346 !--- Pond growing ---! 347 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 348 ! 349 !--- Lid melting ---! 350 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 351 ! 352 !--- mass flux ---! 353 ! MV I would recommend to remove that 354 ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 355 356 IF( zdv_mlt > 0._wp ) THEN 357 zfac = zdv_mlt * rhow * r1_Dt_ice ! melt pond mass flux < 0 [kg.m-2.s-1] 358 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 359 ! 360 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 361 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 362 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 363 ENDIF 364 365 !-------------------! 366 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 367 !-------------------! 368 ! 369 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 370 ! 371 !--- Pond contraction (due to refreezing) ---! 372 zvold = v_ip_1d(ji) ! for diag 373 374 IF( ln_pnd_lids ) THEN 375 ! 376 !--- Lid growing and subsequent pond shrinking ---! 377 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 378 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 379 380 ! Lid growing 381 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 382 383 ! Pond shrinking 384 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 385 386 ELSE 387 388 ! Pond shrinking 389 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 390 391 ENDIF 392 393 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_Dt_ice ! shrinking counted as lid diagnostic 394 395 ! 396 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 397 ! v_ip = h_ip * a_ip 398 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 399 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 400 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 401 402 !------------------------------------------------! 403 ! Pond drainage through brine network (flushing) ! 404 !------------------------------------------------! 405 ! height of top of the pond above sea-level 406 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 407 408 ! Calculate the permeability of the ice 409 ! mv- note here that a linear expression for permeability is fine, 410 ! simpler and more consistent with the rest of SI3 code 411 DO jk = 1, nlay_i 412 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 413 END DO 414 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 415 416 ! Do the drainage using Darcy's law 417 zdv_flush = - zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 418 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 419 ! zdv_flush = 0._wp ! MV remove pond drainage for now 420 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 421 422 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_Dt_ice ! shrinking counted as lid diagnostic 423 424 ! MV --- why pond drainage does not give back water into freshwater flux ? 425 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 426 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 427 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 428 429 !--- Corrections and lid thickness ---! 430 IF( ln_pnd_lids ) THEN 431 !--- retrieve lid thickness from volume ---! 432 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 433 ELSE ; h_il_1d(ji) = 0._wp 434 ENDIF 435 !--- remove ponds if lids are much larger than ponds ---! 436 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 437 a_ip_1d(ji) = 0._wp 438 h_ip_1d(ji) = 0._wp 439 h_il_1d(ji) = 0._wp 440 ENDIF 441 ENDIF 442 ! 443 ENDIF 444 445 END DO ! ji 446 447 !----------------------------------------------------------------------------------------------- 448 ! Retrieve 2D arrays 449 !----------------------------------------------------------------------------------------------- 450 451 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 452 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 453 454 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 455 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 456 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 457 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 458 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 459 460 END DO ! ji 461 462 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 465 466 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 467 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 469 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 470 470 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 471 471 ! 472 ENDIF 473 472 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 473 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 474 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 476 ! 474 477 END SUBROUTINE pnd_LEV 475 478 … … 507 510 !! 508 511 !!------------------------------------------------------------------- 512 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 513 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 514 zvp_min = 1.e-4_wp ! minimum pond volume (m) 515 509 516 510 517 ! local variables … … 536 543 ! a_ip -> apond 537 544 ! a_ip_frac -> apnd 545 546 CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 538 547 539 548 !--------------------------------------------------------------- … … 628 637 DO_2D( 1, 1, 1, 1 ) 629 638 630 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) ) THEN639 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 631 640 632 641 !-------------------------- … … 774 783 775 784 END DO ! jl 785 776 786 777 787 END SUBROUTINE pnd_TOPO … … 848 858 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 849 859 850 INTEGER :: ji, jj, jk, jl ! loop indices 860 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 861 zvp_min = 1.e-4_wp ! minimum pond volume (m) 862 863 INTEGER :: ji, jj, jk, jl ! loop indices 851 864 852 865 a_ip(:,:,:) = 0._wp … … 855 868 DO_2D( 1, 1, 1, 1 ) 856 869 857 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN870 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 858 871 859 872 !------------------------------------------------------------------- -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceupdate.F90
r13643 r13957 94 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 95 95 !!--------------------------------------------------------------------- 96 IF( ln_timing ) CALL timing_start('ice _update')96 IF( ln_timing ) CALL timing_start('iceupdate') 97 97 98 98 IF( kt == nit000 .AND. lwp ) THEN … … 154 154 ! ice-ocean mass flux 155 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 157 157 158 158 ! snw-ocean mass flux … … 160 160 161 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 164 164 165 165 ! Salt flux at the ocean surface … … 172 172 snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step 173 173 ! ! new mass per unit area 174 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )174 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) ) 175 175 ! ! time evolution of snow+ice mass 176 176 snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice … … 286 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice _update')! timing288 IF( ln_timing ) CALL timing_stop ('iceupdate') ! timing 289 289 ! 290 290 END SUBROUTINE ice_update_flx … … 324 324 REAL(wp) :: zflagi ! - - 325 325 !!--------------------------------------------------------------------- 326 IF( ln_timing ) CALL timing_start('ice_update _tau')326 IF( ln_timing ) CALL timing_start('ice_update') 327 327 328 328 IF( kt == nit000 .AND. lwp ) THEN … … 376 376 CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp ) ! lateral boundary condition 377 377 ! 378 IF( ln_timing ) CALL timing_stop('ice_update _tau')378 IF( ln_timing ) CALL timing_stop('ice_update') 379 379 ! 380 380 END SUBROUTINE ice_update_tau -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90
r13911 r13957 236 236 z1_zhmax = 1._wp / hi_max(jpl) 237 237 WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 238 h_i (:,:,jpl) = zhmax238 h_i (:,:,jpl) = zhmax 239 239 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 240 240 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) … … 252 252 ELSEWHERE( h_il(:,:,:) >= zhl_max ) ; a_ip_eff(:,:,:) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow 253 253 ELSEWHERE ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * & ! lid is in between. Expose part of the pond 254 & ( h_il(:,:,:) - zhl_min) / ( zhl_max - zhl_min )254 & ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 255 255 END WHERE 256 256 ! … … 534 534 DO_2D( 1, 1, 1, 1 ) 535 535 ! update exchanges with ocean 536 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice 537 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice 538 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice 536 sfx_res(ji,jj) = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice 537 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice 538 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice 539 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 539 540 ! 540 541 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) … … 551 552 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 552 553 v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 554 h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 555 h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 553 556 ! 554 557 END_2D 555 558 ! 556 !-----------------------------------------------------------------557 ! zap small ponds558 !-----------------------------------------------------------------559 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN560 DO_2D( 1, 1, 1, 1 )561 IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN562 a_ip(ji,jj,jl) = 0._wp563 a_ip_frac(ji,jj,jl) = 0._wp564 v_ip(ji,jj,jl) = 0._wp565 h_ip(ji,jj,jl) = 0._wp566 v_il(ji,jj,jl) = 0._wp567 h_il(ji,jj,jl) = 0._wp568 ENDIF569 END_2D570 ENDIF571 559 END DO 572 560 … … 650 638 psv_i (ji,jj,jl) = 0._wp 651 639 ENDIF 640 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 641 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 642 pv_il (ji,jj,jl) = 0._wp 643 ENDIF 644 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 645 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 646 pv_ip (ji,jj,jl) = 0._wp 647 ENDIF 652 648 END_2D 653 649 ! … … 658 654 WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp 659 655 WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp 660 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)661 WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok662 656 ! 663 657 END SUBROUTINE ice_var_zapneg
Note: See TracChangeset
for help on using the changeset viewer.