Changeset 13809
- Timestamp:
- 2020-11-18T11:04:09+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90
r13284 r13809 34 34 INTEGER :: nice_pnd ! choice of the type of pond scheme 35 35 ! ! associated indices: 36 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 38 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 36 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 38 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 39 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 39 40 40 41 !! * Substitutions … … 60 61 ! 61 62 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 63 ! 64 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 62 65 ! 63 66 END SELECT … … 153 156 !! Holland et al (J. Clim, 2012) 154 157 !!------------------------------------------------------------------- 158 155 159 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 156 160 !! … … 170 174 !! 171 175 INTEGER :: ji, jk ! loop indices 176 172 177 !!------------------------------------------------------------------- 178 173 179 z1_rhow = 1._wp / rhow 174 180 z1_aspect = 1._wp / zaspect … … 196 202 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 197 203 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 198 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?204 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! >0, max for roundoff errors? 199 205 ! 200 206 !--- overflow ---! … … 208 214 ! h_ip_max = 0.5 * h_i 209 215 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 210 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 216 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 211 217 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 212 218 … … 216 222 !--- Lid melting ---! 217 223 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 224 218 225 ! 219 226 !--- mass flux ---! 220 227 IF( zdv_mlt > 0._wp ) THEN 228 ! MV add comment on what that mass flux means 229 ! water removed from fw flux due to melt pond growth 221 230 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 222 231 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 223 232 ! 233 ! MV 234 ! why surface melt and snow fluxes must be adjusted is not clear 235 ! sounds like things are counted twice 236 ! 224 237 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 225 238 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) … … 264 277 265 278 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 279 ! MV the expression for zsbr has wrong sign!!!! 280 ! MV note - the sign-corrected expression is inconsistent 281 ! with the rest of the SI3 code which assumes linear liquidus 282 ! best expression (most consistent for SI3 should be) 283 ! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 266 284 DO jk = 1, nlay_i 267 285 zsbr = - 1.2_wp & … … 269 287 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 270 288 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 271 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 289 290 ztmp(jk) = sz_i_1d(ji,jk) / zsbr ! MV -> this is brine fraction and as it reads is always negative 272 291 END DO 273 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 292 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! MV at present since ztmp is negative, this is always zero 274 293 275 294 ! Do the drainage using Darcy's law … … 301 320 ! 302 321 END SUBROUTINE pnd_LEV 303 322 323 SUBROUTINE pnd_TOPO (aice, aicen, & 324 vice, vicen, & 325 vsnon, & 326 ticen, salin, & 327 a_ip_frac, h_ip, & 328 Tsfc ) 329 330 !!------------------------------------------------------------------- 331 !! *** ROUTINE pnd_TOPO *** 332 !! 333 !! ** Purpose : Compute melt pond evolution 334 !! 335 !! ** Purpose : Compute melt pond evolution based on the ice 336 !! topography as inferred from the ice thickness 337 !! distribution. 338 !! 339 !! ** Method : This code is initially based on Flocco and Feltham 340 !! (2007) and Flocco et al. (2010). More to come... 341 !! 342 !! ** Tunable parameters : 343 !! 344 !! ** Note : 345 !! 346 !! ** References 347 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 348 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 349 !! 10.1029/2006JC003836. 350 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 351 !! a physically based melt pond scheme into the sea ice component of a 352 !! climate model. J. Geophys. Res. 115, C08012, 353 !! doi: 10.1029/2009JC005568. 354 !! 355 !!------------------------------------------------------------------- 356 357 !js 190423: the lid on melt ponds appears only in the analog subroutine of CICE 5.1.2 358 359 REAL (wp), DIMENSION (jpi,jpj), & 360 INTENT(IN) :: & 361 aice, & ! total ice area fraction 362 vice ! total ice volume (m) 363 364 REAL (wp), DIMENSION (jpi,jpj,jpl), & 365 INTENT(IN) :: & 366 aicen, & ! ice area fraction, per category 367 vsnon, & ! snow volume, per category (m) 368 vicen ! ice volume, per category (m) 369 370 REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 371 INTENT(IN) :: & 372 ticen, & ! ice temperature per category (K) 373 salin 374 375 REAL (wp), DIMENSION (jpi,jpj,jpl), & 376 INTENT(INOUT) :: & 377 a_ip_frac , & ! pond area fraction of ice, per ice category 378 h_ip ! pond depth, per ice category (m) 379 380 REAL (wp), DIMENSION (jpi,jpj,jpl), & 381 INTENT(IN) :: & 382 Tsfc ! snow/sea ice surface temperature 383 384 ! local variables 385 REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 386 zTsfcn, & ! ice/snow surface temperature (C) 387 zvolpn ! pond volume per unit area, per category (m) 388 389 REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 390 zrfrac, & ! fraction of available meltwater retained for melt ponding 391 zapondn,& ! pond area fraction, per category 392 zhpondn ! pond depth, per category (m) 393 394 REAL (wp), DIMENSION (jpi,jpj) :: & 395 zvolp, & ! total volume of pond, per unit area of pond (m) 396 zwfx_tmp ! temporary array for melt water 397 398 REAL (wp) :: & 399 zhi, & ! ice thickness (m) 400 zTavg, & ! mean surface temperature across categories (C) 401 z1_rhow, & ! inverse freshwater density 402 zdTs, & ! temperature difference for freeze-up (C) 403 zvpold, & ! dummy pond volume 404 zdvn ! change in melt pond volume for fresh water budget 405 406 INTEGER, DIMENSION (jpi*jpj) :: & 407 indxi, indxj ! compressed indices for cells with ice melting 408 409 INTEGER :: ij,icells,ji,jj,jl ! loop indices 410 411 REAL (wp), PARAMETER :: & 412 zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 413 zTp = -0.15_wp, & ! pond freezing temperature (C) 414 zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 415 zrexp = 0.01_wp, & ! constant melt pond freeze-up rate 416 z01 = 0.01_wp, & 417 z25 = 0.25_wp, & 418 z5 = 0.5_wp 419 420 z1_rhow = 1. / rhow 421 422 !--------------------------------------------------------------- 423 ! Initialization 424 !--------------------------------------------------------------- 425 zhpondn (:,:,:) = 0._wp 426 zapondn (:,:,:) = 0._wp 427 zvolpn (:,:,:) = 0._wp 428 429 zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0 ! Convert in Celsius 430 431 IF ( ln_pnd_fw ) THEN 432 v_ip_b(:,:,:) = v_ip(:,:,:) 433 ELSE 434 v_ip_b(:,:,:) = 0._wp 435 ENDIF 436 437 !------------------------------------------------------------------ 438 ! Available melt water for melt ponding and corresponding fraction 439 !------------------------------------------------------------------ 440 !js 03/05/19 unset restriction on sign of wfx_pnd_in; mask values close to zero for future division 441 !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp ) ! available meltwater for melt ponding 442 !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) * MAX( 0._wp, SIGN( 1._wp, wfx_sum(:,:) + wfx_snw_sum(:,:) - epsi10 ) ) 443 !wfx_pnd_in(:,:) = (wfx_sum(:,:) + wfx_snw_sum(:,:)) * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) 444 445 ! MV 446 ! NB: wfx_pnd_in can be slightly negative for very small values (why?) 447 ! This can in some occasions give negative 448 ! v_ip in the first category, which then gives crazy pond 449 ! fractions and crashes the code as soon as the melt-pond 450 ! radiative coupling is activated 451 ! if we understand and remove why wfx_sum or wfx_snw could be 452 ! negative, then, we can remove the MAX 453 ! NB: I now changed to wfx_snw_sum, this may fix the problem. 454 ! We should check 455 456 457 ! OLI 07/2017: when we (MV & OL) first started the inclusion of melt 458 ! ponds in the model, we removed the Holland et al. (2012, see 459 ! CESM scheme above) parameterization. I put it back here, 460 ! because I think it is needed. In summary, the sinks of FW for 461 ! ponds are: 1/ Runoff through cracks/leads => depends on the 462 ! total ice area only 463 ! 2/ overflow, including Lüthje et al. (2006) limitation 464 ! (max a_ip fraction function of h_i). This is in 465 ! fact an other form of runoff that depends on the 466 ! ITD 467 ! 3/ Flushing - losses by permeability 468 ! 4/ Refreezing 469 ! 5/ Removal of ponds on thin ice 470 ! I think 1 is needed because it is different from 2. However, 471 ! test runs should/could be done, to check the sensitivity and 472 ! the real usefulness of that stuff. 473 ! Note : the Holland et al. param was wrongly wired in NEMO3.1 (using 474 ! a_i instead of at_i), which might well explain why I had a too 475 ! weak melt pond cover in my simulations (compared to MODIS, in 476 ! situ obs. and CICE simulations. 477 478 !js 23/04/19: rewired back to a fraction with a_i 479 zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:) 480 zwfx_tmp(:,:) = 0._wp 481 482 ! MV ---> use expression from level-ice ponds, clarify what is needed 483 484 !--- Add retained melt water to melt ponds 485 ! v_ip should never be negative, otherwise code crashes 486 ! Rq MV: as far as I saw, UM5 can create very small negative v_ip values 487 ! hence I added the max, which was not required with Prather (1 yr run) 488 ! OLI: Here I use vt_ip, so I don't know if the max is 489 ! required... 490 !zvolp(:,:) = MAX(vt_ip(:,:),0._wp) + zrfrac(:,:) * wfx_pnd_in(:,:) * z1_rhow * rdt_ice ! Total available melt water, to be distributed as melt ponds 491 ! OLI: 07/2017 Bugfix above, removed " * aice(:,:)" 492 !js 19/04/18: change zrfrac to use aicen 493 494 zvolp(:,:) = vt_ip(:,:) 495 496 DO jl = 1, jpl 497 ! Melt water, to be distributed as melt ponds 498 zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl) & 499 * ( dh_i_pnd(:,:,jl)*rhoic + dh_s_pnd(:,:,jl)*rhosn ) & 500 * z1_rhow * a_i(:,:,jl) 501 END DO 502 503 ! MV ---> use expression from level ice melt ponds (dv_mlt) 504 505 !js 03/05/19: we truncate negative values after calculating zvolp, in a 506 ! similar manner to the subroutine lim_mp_cesm. Variation dh_i_pnd and 507 ! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them 508 ! to be negative for some reasons. We keep this behaviour as it is, for 509 ! fluxes conservation reasons. If some dh are positive, then we remove water 510 ! indirectly from the ponds. 511 zvolp(:,:) = MAX( zvolp(:,:) , 0._wp ) 512 513 514 ! Fresh water flux going into the ponds 515 wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice 516 517 !--- Remove retained meltwater from surface fluxes 518 IF ( ln_pnd_fw ) THEN 519 !wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrfrac(:,:) ) 520 !wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrfrac(:,:) ) 521 522 !js 190419: we change the code to use a_i in zrfrac. To be tested, but 523 !it should be conservative. zwfx_tmp is the flux accumulated in the 524 !ponds. wfx_pnd_in is the total surface melt fluxes. 525 zwfx_tmp(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) & 526 & * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) 527 WHERE ( ABS(zwfx_tmp(:,:)) > epsi10 ) 528 zwfx_tmp(:,:) = wfx_pnd_in(:,:) / zwfx_tmp(:,:) 529 ELSEWHERE 530 zwfx_tmp(:,:) = 0._wp 531 ENDWHERE 532 533 wfx_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_sum(:,:) 534 wfx_snw_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_snw_sum(:,:) 535 ENDIF 536 537 !----------------------------------------------------------------- 538 ! Identify grid cells with ponds 539 !----------------------------------------------------------------- 540 541 icells = 0 542 DO jj = 1, jpj 543 DO ji = 1, jpi 544 IF ( aice(ji,jj) > epsi10 ) THEN 545 zhi = vice(ji,jj) / aice(ji,jj) 546 ELSE 547 zhi = 0._wp 548 END IF 549 550 IF ( aice(ji,jj) > z01 .and. zhi > rn_himin .and. & 551 zvolp(ji,jj) > zmin_volp*aice(ji,jj)) THEN 552 icells = icells + 1 553 indxi(icells) = ji 554 indxj(icells) = jj 555 ELSE ! remove ponds on thin ice, or too small ponds 556 zvolpn(ji,jj,:) = 0._wp 557 zvolp (ji,jj) = 0._wp 558 559 a_ip(ji,jj,:) = 0._wp 560 v_ip(ji,jj,:) = 0._wp 561 a_ip_frac(ji,jj,:) = 0._wp 562 h_ip(ji,jj,:) = 0._wp 563 564 vt_ip(ji,jj) = 0._wp 565 at_ip(ji,jj) = 0._wp 566 567 ! IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean 568 ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zvolp(ji,jj) * rhow * r1_rdtice 569 END IF 570 END DO ! ji 571 END DO ! jj 572 573 574 DO ij = 1, icells 575 576 ji = indxi(ij) 577 jj = indxj(ij) 578 579 !-------------------------------------------------------------- 580 ! Shrink pond due to refreezing 581 !-------------------------------------------------------------- 582 ! OLI 07/2017: Done like for empirical melt pond scheme (CESM). 583 ! Therefore, I chose to put this part of the code before the main 584 ! routines lim_mp_area/depth (contrary to the original code), 585 ! seeing the freeze-up as a global sink of 586 ! freshwater for melt ponds in the whole grid cell. If this was done 587 ! after, I would need to make an additional assumption on the shape of 588 ! melt ponds, which I don't want to do (for the CESM scheme, this 589 ! assumption was on the aspect ratio). So I remove some water due to 590 ! refreezing first (using zTavg instead of zTsfcn in each category) and 591 ! then let the FF07 routines do their job for the fractional areas and 592 ! depths of melt ponds. 593 ! The whole ice lid related stuff from FF07 was thus removed and replaced 594 ! by this. As mentionned below, this should be improved, but is much 595 ! easier to conserve heat and freshwater this way. 596 597 ! Average surface temperature is needed to compute freeze-up at the cell 598 ! scale 599 zTavg = 0._wp 600 DO jl = 1, jpl 601 zTavg = zTavg + zTsfcn(ji,jj,jl)*aicen(ji,jj,jl) 602 END DO 603 zTavg = zTavg / aice(ji,jj) 604 605 ! The freezing temperature for meltponds is assumed slightly below 0C, 606 ! as if meltponds had a little salt in them (hence the use of zTp). 607 ! The salt budget is not 608 ! altered for meltponds, but if it were then an actual pond freezing 609 ! temperature could be computed. 610 611 zdTs = MAX ( zTp - zTavg, 0. ) 612 613 zvpold = zvolp(ji,jj) 614 615 zvolp(ji,jj) = zvolp(ji,jj) * EXP( zrexp * zdTs / zTp ) 616 617 !--- Dump meltwater due to refreezing ( of course this is wrong 618 !--- but this parameterization is too simple ) 619 ! IF ( ln_pnd_fw ) & 620 ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + rhow * ( zvpold - zvolp(ji,jj) ) * r1_rdtice 621 ! ! OLI 07/2017 : Bugfix above, zvpold - zvolp instead of the 622 ! ! opposite, otherwise negative contribution 623 624 !-------------------------------------------------------------- 625 ! calculate pond area and depth 626 !-------------------------------------------------------------- 627 zdvn = 0._wp 628 CALL lim_mp_area(aice(ji,jj),vice(ji,jj), & 629 aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), & 630 ticen(ji,jj,:,:), salin(ji,jj,:,:), & 631 zvolpn(ji,jj,:), zvolp(ji,jj), & 632 zapondn(ji,jj,:),zhpondn(ji,jj,:), zdvn) 633 ! outputs are 634 ! - zdvn 635 ! - zvolpn 636 ! - zvolp 637 ! - zapondn 638 ! - zhpondn 639 640 ! IF ( ln_pnd_fw ) & 641 ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zdvn * rhow * r1_rdtice ! update flux from ponds to ocean 642 643 !--------------------------------------------------------------- 644 ! Update pond volume and fraction 645 !--------------------------------------------------------------- 646 DO jl = 1, jpl 647 a_ip(ji,jj,jl) = zapondn(ji,jj,jl) 648 v_ip(ji,jj,jl) = zvolpn(ji,jj,jl) 649 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / MAX(aicen(ji,jj,jl), epsi10) & 650 * MAX(0._wp, SIGN(1._wp, aicen(ji,jj,jl) - epsi10)) 651 h_ip (ji,jj,jl) = zhpondn(ji,jj,jl) 652 END DO 653 END DO ! ij 654 655 IF ( ln_pnd_fw ) THEN 656 !js 15/05/19: water going out of the ponds give a positive freshwater 657 ! flux. 658 wfx_pnd_out(:,:) = SUM(MAX(0._wp, v_ip_b(:,:,:) - v_ip(:,:,:)), DIM=3) * rhow * r1_rdtice 659 ELSE 660 wfx_pnd_out(:,:) = 0._wp 661 ENDIF 662 663 END SUBROUTINE pnd_TOPO 664 665 SUBROUTINE lim_mp_area(aice, vice, & 666 aicen, vicen, vsnon, ticen, & 667 salin, zvolpn, zvolp, & 668 zapondn,zhpondn,dvolp) 669 670 !!------------------------------------------------------------------- 671 !! *** ROUTINE lim_mp_area *** 672 !! 673 !! ** Purpose : Given the total volume of meltwater, update 674 !! pond fraction (a_ip) and depth (should be volume) 675 !! 676 !! ** 677 !! 678 !!------------------------------------------------------------------ 679 680 REAL (wp), INTENT(IN) :: & 681 aice, vice 682 683 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 684 aicen, vicen, vsnon 685 686 REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 687 ticen, salin 688 689 REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 690 zvolpn 691 692 REAL (wp), INTENT(INOUT) :: & 693 zvolp, dvolp 694 695 REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 696 zapondn, zhpondn 697 698 INTEGER :: & 699 n, ns, & 700 m_index, & 701 permflag 702 703 REAL (wp), DIMENSION(jpl) :: & 704 hicen, & 705 hsnon, & 706 asnon, & 707 rhos, & ! OLI 07/2017 : for now this is useless, but will be useful with new snow scheme 708 alfan, & 709 betan, & 710 cum_max_vol, & 711 reduced_aicen 712 713 REAL (wp), DIMENSION(0:jpl) :: & 714 cum_max_vol_tmp 715 716 REAL (wp) :: & 717 hpond, & 718 drain, & 719 floe_weight, & 720 pressure_head, & 721 hsl_rel, & 722 deltah, & 723 perm, & 724 msno 725 726 REAL (wp), parameter :: & 727 viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s 728 z0 = 0.0_wp , & 729 c1 = 1.0_wp , & 730 p4 = 0.4_wp , & 731 p6 = 0.6_wp 732 733 !-----------| 734 ! | 735 ! |-----------| 736 !___________|___________|______________________________________sea-level 737 ! | | 738 ! | |---^--------| 739 ! | | | | 740 ! | | | |-----------| |------- 741 ! | | |alfan(n)| | | 742 ! | | | | |--------------| 743 ! | | | | | | 744 !---------------------------v------------------------------------------- 745 ! | | ^ | | | 746 ! | | | | |--------------| 747 ! | | |betan(n)| | | 748 ! | | | |-----------| |------- 749 ! | | | | 750 ! | |---v------- | 751 ! | | 752 ! |-----------| 753 ! | 754 !-----------| 755 756 !------------------------------------------------------------------- 757 ! initialize 758 !------------------------------------------------------------------- 759 760 rhos(:) = rhosn ! OLI 07/2017 : same has above 761 762 DO n = 1, jpl 763 764 zapondn(n) = z0 765 zhpondn(n) = z0 766 767 !---------------------------------------- 768 ! X) compute the effective snow fraction 769 !---------------------------------------- 770 IF (aicen(n) < epsi10) THEN 771 hicen(n) = z0 772 hsnon(n) = z0 773 reduced_aicen(n) = z0 774 asnon(n) = z0 !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 775 ELSE 776 hicen(n) = vicen(n) / aicen(n) 777 hsnon(n) = vsnon(n) / aicen(n) 778 reduced_aicen(n) = c1 ! n=jpl 779 780 !js: initial code in NEMO_DEV 781 !IF (n < jpl) reduced_aicen(n) = aicen(n) & 782 ! * (-0.024_wp*hicen(n) + 0.832_wp) 783 784 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 785 IF (n < jpl) reduced_aicen(n) = aicen(n) & 786 * max(0.2_wp,(-0.024_wp*hicen(n) + 0.832_wp)) 787 788 asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical) 789 ! MV should check whether this makes sense to have the same effective snow fraction in here 790 ! OLI: it probably doesn't 791 END IF 792 793 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 794 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 795 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 796 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 797 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 798 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 799 ! alfan = 60% of the ice volume) in each category lies above the reference line, 800 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 801 802 ! MV: 803 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 804 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 805 806 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 807 808 alfan(n) = 0.6 * hicen(n) 809 betan(n) = 0.4 * hicen(n) 810 811 cum_max_vol(n) = z0 812 cum_max_vol_tmp(n) = z0 813 814 END DO ! jpl 815 816 cum_max_vol_tmp(0) = z0 817 drain = z0 818 dvolp = z0 819 820 !---------------------------------------------------------- 821 ! x) Drain overflow water, update pond fraction and volume 822 !---------------------------------------------------------- 823 824 !-------------------------------------------------------------------------- 825 ! the maximum amount of water that can be contained up to each ice category 826 !-------------------------------------------------------------------------- 827 828 ! MV 829 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 830 ! Then the excess volume cum_max_vol(jl) drains out of the system 831 ! It should be added to wfx_pnd_out 832 ! END MV 833 !js 18/04/19: XXX do something about this flux thing 834 835 DO n = 1, jpl-1 ! last category can not hold any volume 836 837 IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 838 839 ! total volume in level including snow 840 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 841 (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 842 843 ! subtract snow solid volumes from lower categories in current level 844 DO ns = 1, n 845 cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 846 - rhos(ns)/rhow * & ! free air fraction that can be filled by water 847 asnon(ns) * & ! effective areal fraction of snow in that category 848 max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), z0) 849 END DO 850 851 ELSE ! assume higher categories unoccupied 852 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 853 END IF 854 !IF (cum_max_vol_tmp(n) < z0) THEN 855 ! CALL abort_ice('negative melt pond volume') 856 !END IF 857 END DO 858 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 859 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 860 861 !---------------------------------------------------------------- 862 ! is there more meltwater than can be held in the floe? 863 !---------------------------------------------------------------- 864 IF (zvolp >= cum_max_vol(jpl)) THEN 865 drain = zvolp - cum_max_vol(jpl) + epsi10 866 zvolp = zvolp - drain ! update meltwater volume available 867 dvolp = drain ! this is the drained water 868 IF (zvolp < epsi10) THEN 869 dvolp = dvolp + zvolp 870 zvolp = z0 871 END IF 872 END IF 873 874 ! height and area corresponding to the remaining volume 875 876 CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 877 878 DO n=1, m_index 879 !zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 880 ! ! volume instead, no ? 881 zhpondn(n) = max((hpond - alfan(n) + alfan(1)), z0) !js: from CICE 5.1.2 882 zapondn(n) = reduced_aicen(n) 883 ! in practise, pond fraction depends on the empirical snow fraction 884 ! so in turn on ice thickness 885 END DO 886 !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 887 888 !------------------------------------------------------------------------ 889 ! Drainage through brine network (permeability) 890 !------------------------------------------------------------------------ 891 !!! drainage due to ice permeability - Darcy's law 892 893 ! sea water level 894 msno = z0 895 DO n=1,jpl 896 msno = msno + vsnon(n) * rhos(n) 897 END DO 898 floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 899 hsl_rel = floe_weight / rau0 & 900 - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 901 902 deltah = hpond - hsl_rel 903 pressure_head = grav * rau0 * max(deltah, z0) 904 905 ! drain if ice is permeable 906 permflag = 0 907 IF (pressure_head > z0) THEN 908 DO n = 1, jpl-1 909 IF (hicen(n) /= z0) THEN 910 !IF (hicen(n) > z0) THEN !js: from CICE 5.1.2 911 perm = 0._wp ! MV ugly dummy patch 912 CALL lim_mp_perm(ticen(:,n), salin(:,n), perm) 913 IF (perm > z0) permflag = 1 914 915 drain = perm*zapondn(n)*pressure_head*rdt_ice / & 916 (viscosity*hicen(n)) 917 dvolp = dvolp + min(drain, zvolp) 918 zvolp = max(zvolp - drain, z0) 919 IF (zvolp < epsi10) THEN 920 dvolp = dvolp + zvolp 921 zvolp = z0 922 END IF 923 END IF 924 END DO 925 926 ! adjust melt pond dimensions 927 IF (permflag > 0) THEN 928 ! recompute pond depth 929 CALL lim_mp_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 930 DO n=1, m_index 931 zhpondn(n) = hpond - alfan(n) + alfan(1) 932 zapondn(n) = reduced_aicen(n) 933 END DO 934 !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 935 END IF 936 END IF ! pressure_head 937 938 !------------------------------- 939 ! X) remove water from the snow 940 !------------------------------- 941 !------------------------------------------------------------------------ 942 ! total melt pond volume in category does not include snow volume 943 ! snow in melt ponds is not melted 944 !------------------------------------------------------------------------ 945 946 ! Calculate pond volume for lower categories 947 DO n=1,m_index-1 948 zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 949 - (rhos(n)/rhow) * asnon(n) * min(hsnon(n), zhpondn(n)) 950 END DO 951 952 ! Calculate pond volume for highest category = remaining pond volume 953 954 ! The following is completely unclear to Martin at least 955 ! Could we redefine properly and recode in a more readable way ? 956 957 ! m_index = last category with melt pond 958 959 IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 960 961 IF (m_index > 1) THEN 962 IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 963 zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) 964 ELSE 965 zvolpn(m_index) = z0 966 zhpondn(m_index) = z0 967 zapondn(m_index) = z0 968 ! If remaining pond volume is negative reduce pond volume of 969 ! lower category 970 IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 971 zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp 972 END IF 973 END IF 974 975 DO n=1,m_index 976 IF (zapondn(n) > epsi10) THEN 977 zhpondn(n) = zvolpn(n) / zapondn(n) 978 ELSE 979 dvolp = dvolp + zvolpn(n) 980 zhpondn(n) = z0 981 zvolpn(n) = z0 982 zapondn(n) = z0 983 END IF 984 END DO 985 DO n = m_index+1, jpl 986 zhpondn(n) = z0 987 zapondn(n) = z0 988 zvolpn (n) = z0 989 END DO 990 991 END SUBROUTINE lim_mp_area 992 993 994 SUBROUTINE lim_mp_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index) 995 !!------------------------------------------------------------------- 996 !! *** ROUTINE lim_mp_depth *** 997 !! 998 !! ** Purpose : Compute melt pond depth 999 !!------------------------------------------------------------------- 1000 1001 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1002 aicen, & 1003 asnon, & 1004 hsnon, & 1005 rhos, & 1006 alfan, & 1007 cum_max_vol 1008 1009 REAL (wp), INTENT(IN) :: & 1010 zvolp 1011 1012 REAL (wp), INTENT(OUT) :: & 1013 hpond 1014 1015 INTEGER, INTENT(OUT) :: & 1016 m_index 1017 1018 INTEGER :: n, ns 1019 1020 REAL (wp), DIMENSION(0:jpl+1) :: & 1021 hitl, & 1022 aicetl 1023 1024 REAL (wp) :: & 1025 rem_vol, & 1026 area, & 1027 vol, & 1028 tmp, & 1029 z0 = 0.0_wp 1030 1031 !---------------------------------------------------------------- 1032 ! hpond is zero if zvolp is zero - have we fully drained? 1033 !---------------------------------------------------------------- 1034 1035 IF (zvolp < epsi10) THEN 1036 hpond = z0 1037 m_index = 0 1038 ELSE 1039 1040 !---------------------------------------------------------------- 1041 ! Calculate the category where water fills up to 1042 !---------------------------------------------------------------- 1043 1044 !----------| 1045 ! | 1046 ! | 1047 ! |----------| -- -- 1048 !__________|__________|_________________________________________ ^ 1049 ! | | rem_vol ^ | Semi-filled 1050 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1051 ! | | | | 1052 ! | | | |hpond 1053 ! | | |----------| | |------- 1054 ! | | | | | | 1055 ! | | | |---v-----| 1056 ! | | m_index | | | 1057 !------------------------------------------------------------- 1058 1059 m_index = 0 ! 1:m_index categories have water in them 1060 DO n = 1, jpl 1061 IF (zvolp <= cum_max_vol(n)) THEN 1062 m_index = n 1063 IF (n == 1) THEN 1064 rem_vol = zvolp 1065 ELSE 1066 rem_vol = zvolp - cum_max_vol(n-1) 1067 END IF 1068 exit ! to break out of the loop 1069 END IF 1070 END DO 1071 m_index = min(jpl-1, m_index) 1072 1073 !---------------------------------------------------------------- 1074 ! semi-filled layer may have m_index different snow in it 1075 !---------------------------------------------------------------- 1076 1077 !----------------------------------------------------------- ^ 1078 ! | alfan(m_index+1) 1079 ! | 1080 !hitl(3)--> |----------| | 1081 !hitl(2)--> |------------| * * * * *| | 1082 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1083 !hitl(0)-->------------------------------------------------- | ^ 1084 ! various snow from lower categories | |alfa(m_index) 1085 1086 ! hitl - heights of the snow layers from thinner and current categories 1087 ! aicetl - area of each snow depth in this layer 1088 1089 hitl(:) = z0 1090 aicetl(:) = z0 1091 DO n = 1, m_index 1092 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1093 alfan(m_index+1) - alfan(m_index)), z0) 1094 aicetl(n) = asnon(n) 1095 1096 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1097 END DO 1098 1099 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1100 aicetl(m_index+1) = z0 1101 1102 !---------------------------------------------------------------- 1103 ! reorder array according to hitl 1104 ! snow heights not necessarily in height order 1105 !---------------------------------------------------------------- 1106 1107 DO ns = 1, m_index+1 1108 DO n = 0, m_index - ns + 1 1109 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1110 tmp = hitl(n) 1111 hitl(n) = hitl(n+1) 1112 hitl(n+1) = tmp 1113 tmp = aicetl(n) 1114 aicetl(n) = aicetl(n+1) 1115 aicetl(n+1) = tmp 1116 END IF 1117 END DO 1118 END DO 1119 1120 !---------------------------------------------------------------- 1121 ! divide semi-filled layer into set of sublayers each vertically homogenous 1122 !---------------------------------------------------------------- 1123 1124 !hitl(3)---------------------------------------------------------------- 1125 ! | * * * * * * * * 1126 ! |* * * * * * * * * 1127 !hitl(2)---------------------------------------------------------------- 1128 ! | * * * * * * * * | * * * * * * * * 1129 ! |* * * * * * * * * |* * * * * * * * * 1130 !hitl(1)---------------------------------------------------------------- 1131 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1132 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1133 !hitl(0)---------------------------------------------------------------- 1134 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1135 1136 ! move up over layers incrementing volume 1137 DO n = 1, m_index+1 1138 1139 area = sum(aicetl(:)) - & ! total area of sub-layer 1140 (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1141 1142 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1143 1144 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1145 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1146 1147 exit 1148 ELSE ! still in sub-layer below the sub-layer with the depth 1149 rem_vol = rem_vol - vol 1150 END IF 1151 1152 END DO 1153 1154 END IF 1155 1156 END SUBROUTINE lim_mp_depth 1157 1158 1159 SUBROUTINE lim_mp_perm(ticen, salin, perm) 1160 !!------------------------------------------------------------------- 1161 !! *** ROUTINE lim_mp_perm *** 1162 !! 1163 !! ** Purpose : Determine the liquid fraction of brine in the ice 1164 !! and its permeability 1165 !!------------------------------------------------------------------- 1166 1167 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1168 ticen, & ! internal ice temperature (K) 1169 salin ! salinity (ppt) !js: ppt according to cice 1170 1171 REAL (wp), INTENT(OUT) :: & 1172 perm ! permeability 1173 1174 REAL (wp) :: & 1175 Sbr ! brine salinity 1176 1177 REAL (wp), DIMENSION(nlay_i) :: & 1178 Tin, & ! ice temperature 1179 phi ! liquid fraction 1180 1181 INTEGER :: k 1182 1183 !----------------------------------------------------------------- 1184 ! Compute ice temperatures from enthalpies using quadratic formula 1185 !----------------------------------------------------------------- 1186 1187 DO k = 1,nlay_i 1188 Tin(k) = ticen(k) - rt0 !js: from K to degC 1189 END DO 1190 1191 !----------------------------------------------------------------- 1192 ! brine salinity and liquid fraction 1193 !----------------------------------------------------------------- 1194 1195 IF (maxval(Tin) <= -2.0_wp) THEN 1196 1197 ! Assur 1958 1198 DO k = 1,nlay_i 1199 Sbr = - 1.2_wp & 1200 - 21.8_wp * Tin(k) & 1201 - 0.919_wp * Tin(k)**2 & 1202 - 0.01878_wp * Tin(k)**3 1203 phi(k) = salin(k) / Sbr ! liquid fraction 1204 1205 !js: Sbr must not be zero 1206 END DO ! k 1207 1208 ELSE 1209 1210 ! Notz 2005 thesis eq. 3.2 1211 !js: "interstitial brine Sbr (in ppt) as a function of temperature T (in degC)". 1212 DO k = 1,nlay_i 1213 Sbr = - 17.6_wp * Tin(k) & 1214 - 0.389_wp * Tin(k)**2 & 1215 - 0.00362_wp* Tin(k)**3 1216 phi(k) = salin(k) / Sbr ! liquid fraction 1217 1218 !js: Sbr must not be zero 1219 END DO 1220 1221 END IF 1222 1223 !----------------------------------------------------------------- 1224 ! permeability 1225 !----------------------------------------------------------------- 1226 1227 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1228 1229 END SUBROUTINE lim_mp_perm 1230 1231 1232 !---------------------------------------------------------------------------------------------------------------------- 304 1233 305 1234 SUBROUTINE ice_thd_pnd_init … … 319 1248 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 320 1249 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1250 & ln_pnd_TOPO , & 321 1251 & ln_pnd_lids, ln_pnd_alb 322 1252 !!------------------------------------------------------------------- … … 336 1266 WRITE(numout,*) ' Namelist namicethd_pnd:' 337 1267 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 1268 WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO 338 1269 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 339 1270 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min
Note: See TracChangeset
for help on using the changeset viewer.