Changeset 13891
- Timestamp:
- 2020-11-26T17:02:26+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
r13879 r13891 45 45 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 46 46 zvp_min = 1.e-4_wp ! minimum pond volume (m) 47 ! 48 ! Pond volume per area budget diags 49 ! 50 ! The idea of diags is the volume of ponds per grid cell area is 51 ! 52 ! dV/dt = mlt + drn + lid + rnf 53 ! mlt = input from surface melting 54 ! drn = drainage through brine network 55 ! lid = lid growth & melt 56 ! rnf = runoff (water directly removed out of surface melting + overflow) 57 ! 58 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 59 ! 60 ! In level mode 61 62 REAL(wp), DIMENSION(jpi,jpj) :: ! pond volume budget diagnostics 63 diag_dvpn_mlt, & !! meltwater pond volume input [m/day] 64 diag_dvpn_drn, & !! pond volume lost by drainage [m/day] 65 diag_dvpn_lid, & !! pond volume lost by refreezing [m/day] 66 diag_dvpn_rnf !! meltwater pond lost to runoff [m/day] 67 68 REAL(wp), DIMENSION(jpij) :: ! pond volume budget diagnostics (1d) 69 diag_dvpn_mlt_1d, & !! meltwater pond volume input [m/day] 70 diag_dvpn_drn_1d, & !! pond volume lost by drainage [m/day] 71 diag_dvpn_lid_1d, & !! pond volume lost by refreezing [m/day] 72 diag_dvpn_rnf_1d !! meltwater pond lost to runoff [m/day] 73 47 74 48 75 !! * Substitutions … … 76 103 !!------------------------------------------------------------------- 77 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 109 78 110 SELECT CASE ( nice_pnd ) 79 111 ! … … 82 114 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 83 115 ! 84 CASE (np_pndTOPO) ; CALL pnd_TOPO !&!== Topographic melt ponds ==!116 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 85 117 ! 86 118 END SELECT 87 119 ! 120 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 125 88 126 END SUBROUTINE ice_thd_pnd 127 89 128 90 129 … … 191 230 REAL(wp) :: zfac, zdum ! temporary arrays 192 231 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 232 REAL(wp) :: z1_rdtice ! inverse time step 233 REAL(wp) :: zvold ! 193 234 !! 194 235 INTEGER :: ji, jj, jk, jl ! loop indices … … 199 240 z1_aspect = 1._wp / zaspect 200 241 z1_Tp = 1._wp / zTp 242 z1_rdtice = 1._wp / rdt_ice 201 243 202 244 !----------------------------------------------------------------------------------------------- … … 221 263 IF( npti > 0 ) THEN 222 264 223 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 224 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 225 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 265 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 266 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti) , wfx_sum ) 267 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti) , wfx_pnd ) 268 269 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 271 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 272 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 226 273 227 274 DO jl = 1, jpl … … 268 315 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 269 316 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 317 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 270 322 ! 271 323 !--- overflow ---! 272 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 324 ! 325 ! 1) area driven overflow 326 ! 327 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 273 328 ! a_ip_max = zfr_mlt * a_i 274 329 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 275 330 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 276 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 277 331 zvold = zdv_mlt 332 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 333 334 ! 335 ! 2) depth driven overflow 336 ! 278 337 ! If pond depth exceeds half the ice thickness then reduce the pond volume 279 338 ! h_ip_max = 0.5 * h_i 280 339 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 281 340 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 282 342 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 283 345 284 346 !--- Pond growing ---! … … 308 370 ! 309 371 !--- Pond contraction (due to refreezing) ---! 372 zvold = v_ip_1d(ji) ! for diag 373 310 374 IF( ln_pnd_lids ) THEN 311 375 ! … … 319 383 ! Pond shrinking 320 384 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 321 385 322 386 ELSE 387 323 388 ! Pond shrinking 324 389 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 390 325 391 ENDIF 392 393 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * z1_rdtice ! shrinking counted as lid diagnostic 394 326 395 ! 327 396 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac … … 331 400 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 332 401 333 !--------------- !334 ! Pond flushing!335 !--------------- !402 !------------------------------------------------! 403 ! Pond drainage through brine network (flushing) ! 404 !------------------------------------------------! 336 405 ! height of top of the pond above sea-level 337 406 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 … … 352 421 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 353 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 354 zdv_flush = 0._wp ! MV remove pond drainage for now423 ! zdv_flush = 0._wp ! MV remove pond drainage for now 355 424 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 425 426 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * z1_rdtice ! shrinking counted as lid diagnostic 356 427 357 ! MV --- why pond drainage does not give back water into freshwater flux ? 358 428 ! MV --- why pond drainage does not give back water into freshwater flux ? 359 429 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 360 430 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 … … 385 455 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 386 456 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 457 387 458 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 388 459 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) … … 399 470 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 400 471 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 401 472 473 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 474 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 476 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 477 402 478 ! 403 479 ENDIF … … 438 514 ! local variables 439 515 REAL(wp) :: & 440 zdHui, & ! change in thickness of ice lid (m) 441 zomega, & ! conduction 442 zdTice, & ! temperature difference across ice lid (C) 443 zdvice, & ! change in ice volume (m) 444 zTavg, & ! mean surface temperature across categories (C) 445 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 446 zTp, & ! pond freezing temperature (C) 447 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 448 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 449 z1_rhow, & ! inverse water density 450 zpond , & ! dummy variable 451 zdum ! dummy variable 452 453 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 454 zvolp_res 455 516 zdHui, & ! change in thickness of ice lid (m) 517 zomega, & ! conduction 518 zdTice, & ! temperature difference across ice lid (C) 519 zdvice, & ! change in ice volume (m) 520 zTavg, & ! mean surface temperature across categories (C) 521 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 522 zTp, & ! pond freezing temperature (C) 523 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 524 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 525 z1_rhow, & ! inverse water density 526 z1_rdtice, & ! inverse time step 527 zv_pnd , & ! volume of meltwater contributing to ponds 528 zv_mlt ! total amount of meltwater produced 529 530 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 531 zvolp_res !! remaining melt pond water available after drainage 532 456 533 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 457 534 … … 473 550 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 474 551 z1_rhow = 1._wp / rhow 552 z1_rdtice = 1._wp / rdt_ice 475 553 476 554 ! Set required ice variables (hard-coded here for now) … … 491 569 ! Change 2D to 1D 492 570 !--------------------------------------------------------------- 571 572 ! use what we have in iceitd.F90 (incremental remapping) 493 573 494 574 !-------------------------------------------------------------- 495 ! Collect total available pond water 575 ! Collect total available pond water volume 496 576 !-------------------------------------------------------------- 497 577 ! Assuming that meltwater (+rain in principle) runsoff the surface 578 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 579 ! I cite her words, they are very talkative 580 ! "grid cells with very little ice cover (and hence more open water area) 581 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 582 ! "This results in the same runoff fraction r for each ice category within a grid cell" 583 498 584 zvolp(:,:) = 0._wp 499 585 … … 505 591 506 592 !--- Available meltwater for melt ponding ---! 507 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) ! = ( 1 - r ) = fraction of melt water that is not flushed 508 zdum = -( dh_i_sum_2d(ji,jj,jl)*rhoi + dh_s_mlt_2d(ji,jj,jl)*rhos ) * z1_rhow * a_i(ji,jj,jl) 509 zpond = zfr_mlt * zdum ! check 593 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 510 602 511 603 !--- Create possible new ponds … … 513 605 IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 514 606 h_ip(ji,jj,jl) = 0._wp 515 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 607 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) 516 609 ENDIF 517 610 518 !--- Deepen existing ponds before redistribution and drainage(later on) 519 h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl) 520 zvolp(ji,jj) = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 521 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless 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 615 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 616 617 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 618 619 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 522 620 523 621 ENDIF ! a_i … … 526 624 END DO ! jj 527 625 END DO ! ji 528 529 h_ip(:,:,:) = 0._wp ! pond thickness 530 a_ip(:,:,:) = 0._wp ! pond fraction per grid area 531 626 532 627 !-------------------------------------------------------------- 533 628 ! Redistribute and drain water from ponds … … 564 659 565 660 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 566 IF ( zdvice > epsi10 ) then 661 662 IF ( zdvice > epsi10 ) THEN 663 567 664 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 568 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice 665 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 666 ! as it is already counted in surface melt 569 667 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 570 668 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) … … 577 675 ENDIF 578 676 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 677 678 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 679 579 680 ENDIF 580 681 581 682 !---------------------------------------------------------------- 582 683 ! Freeze pre-existing lid … … 598 699 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 599 700 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 701 702 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 703 600 704 ENDIF 601 705 … … 624 728 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 625 729 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 730 731 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 626 732 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 627 733 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice … … 849 955 DO ns = 1, jl 850 956 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 851 - rhos/rhow * & ! free air fraction that can be filled by water957 - rhos/rhow * & ! free air fraction that can be filled by water 852 958 asnon(ns) * & ! effective areal fraction of snow in that category 853 959 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) … … 870 976 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 871 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 872 981 zdvolp(ji,jj) = drain ! this is the drained water 873 982 IF (zvolp(ji,jj) < epsi10) THEN … … 878 987 879 988 ! height and area corresponding to the remaining volume 989 ! routine leaves zvolp unchanged 880 990 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 881 991 … … 917 1027 918 1028 perm = 0._wp ! MV ugly dummy patch 919 !CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof1029 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 920 1030 IF (perm > 0._wp) permflag = 1 921 1031 … … 924 1034 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 925 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 926 1039 IF (zvolp(ji,jj) < epsi10) THEN 927 1040 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) … … 950 1063 ! snow in melt ponds is not melted 951 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 952 1070 953 1071 ! Calculate pond volume for lower categories 954 1072 DO jl = 1,m_index-1 955 1073 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 956 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))1074 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 957 1075 END DO 958 1076
Note: See TracChangeset
for help on using the changeset viewer.
