Changeset 13912
- Timestamp:
- 2020-11-30T11:17:33+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90
r13911 r13912 475 475 END SUBROUTINE pnd_LEV 476 476 477 478 479 SUBROUTINE pnd_TOPO 480 481 !!------------------------------------------------------------------- 482 !! *** ROUTINE pnd_TOPO *** 483 !! 484 !! ** Purpose : Compute melt pond evolution based on the ice 485 !! topography inferred from the ice thickness distribution 486 !! 487 !! ** Method : This code is initially based on Flocco and Feltham 488 !! (2007) and Flocco et al. (2010). 489 !! 490 !! - Calculate available pond water base on surface meltwater 491 !! - Redistribute water as a function of topography, drain water 492 !! - Exchange water with the lid 493 !! 494 !! ** Tunable parameters : 495 !! 496 !! ** Note : 497 !! 498 !! ** References 499 !! 500 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 501 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 502 !! 10.1029/2006JC003836. 503 !! 504 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 505 !! a physically based melt pond scheme into the sea ice component of a 506 !! climate model. J. Geophys. Res. 115, C08012, 507 !! doi: 10.1029/2009JC005568. 508 !! 509 !!------------------------------------------------------------------- 510 511 ! local variables 512 REAL(wp) :: & 513 zdHui, & ! change in thickness of ice lid (m) 514 zomega, & ! conduction 515 zdTice, & ! temperature difference across ice lid (C) 516 zdvice, & ! change in ice volume (m) 517 zTavg, & ! mean surface temperature across categories (C) 518 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 519 zTp, & ! pond freezing temperature (C) 520 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 521 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 522 z1_rhow, & ! inverse water density 523 zv_pnd , & ! volume of meltwater contributing to ponds 524 zv_mlt ! total amount of meltwater produced 525 526 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 527 zvolp_res !! remaining melt pond water available after drainage 528 529 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 530 531 INTEGER :: ji, jj, jk, jl ! loop indices 532 533 INTEGER :: i_test 534 535 ! Note 536 ! equivalent for CICE translation 537 ! a_ip -> apond 538 ! a_ip_frac -> apnd 539 540 !--------------------------------------------------------------- 541 ! Initialise 542 !--------------------------------------------------------------- 543 544 ! Parameters & constants (move to parameters) 545 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 546 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 547 z1_rhow = 1._wp / rhow 548 549 ! Set required ice variables (hard-coded here for now) 550 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 551 552 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 553 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 554 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 555 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 556 557 ! thickness 558 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 559 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 560 END WHERE 561 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 562 563 !--------------------------------------------------------------- 564 ! Change 2D to 1D 565 !--------------------------------------------------------------- 566 ! MV 567 ! a less computing-intensive version would have 2D-1D passage here 568 ! use what we have in iceitd.F90 (incremental remapping) 569 570 !-------------------------------------------------------------- 571 ! Collect total available pond water volume 572 !-------------------------------------------------------------- 573 ! Assuming that meltwater (+rain in principle) runsoff the surface 574 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 575 ! I cite her words, they are very talkative 576 ! "grid cells with very little ice cover (and hence more open water area) 577 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 578 ! "This results in the same runoff fraction r for each ice category within a grid cell" 579 580 zvolp(:,:) = 0._wp 581 582 DO jl = 1, jpl 583 DO_2D( 1, 1, 1, 1 ) 584 585 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 586 587 !--- Available and contributing meltwater for melt ponding ---! 588 zv_mlt = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) & ! available volume of surface melt water per grid area 589 & * z1_rhow * a_i(ji,jj,jl) 590 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 591 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 592 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 593 594 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice ! diags 595 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice 596 597 !--- Create possible new ponds 598 ! if pond does not exist, create new pond over full ice area 599 !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 600 IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 601 a_ip(ji,jj,jl) = a_i(ji,jj,jl) 602 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 603 ENDIF 604 605 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 606 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 607 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 608 609 !--- Total available pond water volume (pre-existing + newly produced)j 610 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 611 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 612 613 ENDIF ! a_i 614 615 END_2D 616 END DO ! ji 617 618 !-------------------------------------------------------------- 619 ! Redistribute and drain water from ponds 620 !-------------------------------------------------------------- 621 CALL ice_thd_pnd_area( zvolp, zvolp_res ) 622 623 !-------------------------------------------------------------- 624 ! Melt pond lid growth and melt 625 !-------------------------------------------------------------- 626 627 IF( ln_pnd_lids ) THEN 628 629 DO_2D( 1, 1, 1, 1 ) 630 631 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 632 633 !-------------------------- 634 ! Pond lid growth and melt 635 !-------------------------- 636 ! Mean surface temperature 637 zTavg = 0._wp 638 DO jl = 1, jpl 639 zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 640 END DO 641 zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 642 643 DO jl = 1, jpl-1 644 645 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 646 647 !---------------------------------------------------------------- 648 ! Lid melting: floating upper ice layer melts in whole or part 649 !---------------------------------------------------------------- 650 ! Use Tsfc for each category 651 IF ( t_su(ji,jj,jl) > zTp ) THEN 652 653 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 654 655 IF ( zdvice > epsi10 ) THEN 656 657 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 658 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 659 ! as it is already counted in surface melt 660 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 661 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 662 663 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 664 ! ice lid melted and category is pond covered 665 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 666 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 667 v_il(ji,jj,jl) = 0._wp 668 ENDIF 669 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 670 671 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 672 673 ENDIF 674 675 !---------------------------------------------------------------- 676 ! Freeze pre-existing lid 677 !---------------------------------------------------------------- 678 679 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 680 681 ! differential growth of base of surface floating ice layer 682 zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 683 zomega = rcnd_i * zdTice / zrhoi_L 684 zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 685 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 686 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 687 688 IF ( zdvice > epsi10 ) THEN 689 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 690 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 691 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 692 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 693 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 694 695 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 696 697 ENDIF 698 699 ENDIF ! Tsfcn(i,j,n) 700 701 !---------------------------------------------------------------- 702 ! Freeze new lids 703 !---------------------------------------------------------------- 704 ! upper ice layer begins to form 705 ! note: albedo does not change 706 707 ELSE ! v_il < epsi10 708 709 ! thickness of newly formed ice 710 ! the surface temperature of a meltpond is the same as that 711 ! of the ice underneath (0C), and the thermodynamic surface 712 ! flux is the same 713 714 !!! we need net surface energy flux, excluding conduction 715 !!! fsurf is summed over categories in CICE 716 !!! we have the category-dependent flux, let us use it ? 717 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 718 zdHui = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) 719 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 720 IF ( zdvice > epsi10 ) THEN 721 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 722 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 723 724 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 725 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 726 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 727 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 728 ENDIF 729 730 ENDIF ! v_il 731 732 END DO ! jl 733 734 ELSE ! remove ponds on thin ice 735 736 v_ip(ji,jj,:) = 0._wp 737 v_il(ji,jj,:) = 0._wp 738 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 739 ! zvolp(ji,jj) = 0._wp 740 741 ENDIF 742 743 END_2D 744 745 ENDIF ! ln_pnd_lids 746 747 !--------------------------------------------------------------- 748 ! Clean-up variables (probably duplicates what icevar would do) 749 !--------------------------------------------------------------- 750 ! MV comment 751 ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 752 ! icevar should recompute all other variables (if needed at all) 753 754 DO jl = 1, jpl 755 756 DO_2D( 1, 1, 1, 1 ) 757 758 ! ! zap lids on small ponds 759 ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 760 ! .AND. v_il(ji,jj,jl) > epsi10) THEN 761 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 762 ! ENDIF 763 764 ! recalculate equivalent pond variables 765 IF ( a_ip(ji,jj,jl) > epsi10) THEN 766 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 767 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 768 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 769 ENDIF 770 ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 771 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 772 ! ENDIF 773 774 END_2D 775 776 END DO ! jl 777 778 END SUBROUTINE pnd_TOPO 779 780 781 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 782 783 !!------------------------------------------------------------------- 784 !! *** ROUTINE ice_thd_pnd_area *** 785 !! 786 !! ** Purpose : Given the total volume of available pond water, 787 !! redistribute and drain water 788 !! 789 !! ** Method 790 !! 791 !-----------| 792 ! | 793 ! |-----------| 794 !___________|___________|______________________________________sea-level 795 ! | | 796 ! | |---^--------| 797 ! | | | | 798 ! | | | |-----------| |------- 799 ! | | | alfan | | | 800 ! | | | | |--------------| 801 ! | | | | | | 802 !---------------------------v------------------------------------------- 803 ! | | ^ | | | 804 ! | | | | |--------------| 805 ! | | | betan | | | 806 ! | | | |-----------| |------- 807 ! | | | | 808 ! | |---v------- | 809 ! | | 810 ! |-----------| 811 ! | 812 !-----------| 813 ! 814 !! 815 !!------------------------------------------------------------------ 816 817 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 818 zvolp, & ! total available pond water 819 zdvolp ! remaining meltwater after redistribution 820 821 INTEGER :: & 822 ns, & 823 m_index, & 824 permflag 825 826 REAL (wp), DIMENSION(jpl) :: & 827 hicen, & 828 hsnon, & 829 asnon, & 830 alfan, & 831 betan, & 832 cum_max_vol, & 833 reduced_aicen 834 835 REAL (wp), DIMENSION(0:jpl) :: & 836 cum_max_vol_tmp 837 838 REAL (wp) :: & 839 hpond, & 840 drain, & 841 floe_weight, & 842 pressure_head, & 843 hsl_rel, & 844 deltah, & 845 perm, & 846 msno 847 848 REAL (wp), parameter :: & 849 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 850 851 INTEGER :: ji, jj, jk, jl ! loop indices 852 853 a_ip(:,:,:) = 0._wp 854 h_ip(:,:,:) = 0._wp 855 856 DO_2D( 1, 1, 1, 1 ) 857 858 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 859 860 !------------------------------------------------------------------- 861 ! initialize 862 !------------------------------------------------------------------- 863 864 DO jl = 1, jpl 865 866 !---------------------------------------- 867 ! compute the effective snow fraction 868 !---------------------------------------- 869 870 IF (a_i(ji,jj,jl) < epsi10) THEN 871 hicen(jl) = 0._wp 872 hsnon(jl) = 0._wp 873 reduced_aicen(jl) = 0._wp 874 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 875 ELSE 876 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 877 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 878 reduced_aicen(jl) = 1._wp ! n=jpl 879 880 !js: initial code in NEMO_DEV 881 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 882 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 883 884 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 885 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 886 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 887 888 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 889 ! MV should check whether this makes sense to have the same effective snow fraction in here 890 ! OLI: it probably doesn't 891 END IF 892 893 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 894 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 895 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 896 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 897 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 898 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 899 ! alfan = 60% of the ice volume) in each category lies above the reference line, 900 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 901 902 ! MV: 903 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 904 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 905 906 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 907 908 alfan(jl) = 0.6 * hicen(jl) 909 betan(jl) = 0.4 * hicen(jl) 910 911 cum_max_vol(jl) = 0._wp 912 cum_max_vol_tmp(jl) = 0._wp 913 914 END DO ! jpl 915 916 cum_max_vol_tmp(0) = 0._wp 917 drain = 0._wp 918 zdvolp(ji,jj) = 0._wp 919 920 !---------------------------------------------------------- 921 ! Drain overflow water, update pond fraction and volume 922 !---------------------------------------------------------- 923 924 !-------------------------------------------------------------------------- 925 ! the maximum amount of water that can be contained up to each ice category 926 !-------------------------------------------------------------------------- 927 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 928 ! Then the excess volume cum_max_vol(jl) drains out of the system 929 ! It should be added to wfx_pnd_out 930 931 DO jl = 1, jpl-1 ! last category can not hold any volume 932 933 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 934 935 ! total volume in level including snow 936 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 937 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 938 939 ! subtract snow solid volumes from lower categories in current level 940 DO ns = 1, jl 941 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 942 - rhos/rhow * & ! free air fraction that can be filled by water 943 asnon(ns) * & ! effective areal fraction of snow in that category 944 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 945 END DO 946 947 ELSE ! assume higher categories unoccupied 948 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 949 END IF 950 !IF (cum_max_vol_tmp(jl) < z0) THEN 951 ! CALL abort_ice('negative melt pond volume') 952 !END IF 953 END DO 954 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 955 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 956 957 !---------------------------------------------------------------- 958 ! is there more meltwater than can be held in the floe? 959 !---------------------------------------------------------------- 960 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 961 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 962 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 963 964 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 965 966 zdvolp(ji,jj) = drain ! this is the drained water 967 IF (zvolp(ji,jj) < epsi10) THEN 968 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 969 zvolp(ji,jj) = 0._wp 970 END IF 971 END IF 972 973 ! height and area corresponding to the remaining volume 974 ! routine leaves zvolp unchanged 975 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 976 977 DO jl = 1, m_index 978 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 979 ! ! volume instead, no ? 980 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 981 a_ip(ji,jj,jl) = reduced_aicen(jl) 982 ! in practise, pond fraction depends on the empirical snow fraction 983 ! so in turn on ice thickness 984 END DO 985 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 986 987 !------------------------------------------------------------------------ 988 ! Drainage through brine network (permeability) 989 !------------------------------------------------------------------------ 990 !!! drainage due to ice permeability - Darcy's law 991 992 ! sea water level 993 msno = 0._wp 994 DO jl = 1 , jpl 995 msno = msno + v_s(ji,jj,jl) * rhos 996 END DO 997 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj) 998 hsl_rel = floe_weight / rho0 & 999 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1000 1001 deltah = hpond - hsl_rel 1002 pressure_head = grav * rho0 * max(deltah, 0._wp) 1003 1004 ! drain if ice is permeable 1005 permflag = 0 1006 1007 IF (pressure_head > 0._wp) THEN 1008 DO jl = 1, jpl-1 1009 IF ( hicen(jl) /= 0._wp ) THEN 1010 1011 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1012 1013 perm = 0._wp ! MV ugly dummy patch 1014 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1015 IF (perm > 0._wp) permflag = 1 1016 1017 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & 1018 (viscosity*hicen(jl)) 1019 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1020 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1021 1022 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1023 1024 IF (zvolp(ji,jj) < epsi10) THEN 1025 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1026 zvolp(ji,jj) = 0._wp 1027 END IF 1028 END IF 1029 END DO 1030 1031 ! adjust melt pond dimensions 1032 IF (permflag > 0) THEN 1033 ! recompute pond depth 1034 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1035 DO jl = 1, m_index 1036 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1037 a_ip(ji,jj,jl) = reduced_aicen(jl) 1038 END DO 1039 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1040 END IF 1041 END IF ! pressure_head 1042 1043 !------------------------------- 1044 ! remove water from the snow 1045 !------------------------------- 1046 !------------------------------------------------------------------------ 1047 ! total melt pond volume in category does not include snow volume 1048 ! snow in melt ponds is not melted 1049 !------------------------------------------------------------------------ 1050 1051 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1052 ! how much, so I did not diagnose it 1053 ! so if there is a problem here, nobody is going to see it... 1054 1055 1056 ! Calculate pond volume for lower categories 1057 DO jl = 1,m_index-1 1058 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1059 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1060 END DO 1061 1062 ! Calculate pond volume for highest category = remaining pond volume 1063 1064 ! The following is completely unclear to Martin at least 1065 ! Could we redefine properly and recode in a more readable way ? 1066 1067 ! m_index = last category with melt pond 1068 1069 IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 1070 1071 IF (m_index > 1) THEN 1072 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1073 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1074 ELSE 1075 v_ip(ji,jj,m_index) = 0._wp 1076 h_ip(ji,jj,m_index) = 0._wp 1077 a_ip(ji,jj,m_index) = 0._wp 1078 ! If remaining pond volume is negative reduce pond volume of 1079 ! lower category 1080 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1081 v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 1082 END IF 1083 END IF 1084 1085 DO jl = 1,m_index 1086 IF (a_ip(ji,jj,jl) > epsi10) THEN 1087 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1088 ELSE 1089 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1090 h_ip(ji,jj,jl) = 0._wp 1091 v_ip(ji,jj,jl) = 0._wp 1092 a_ip(ji,jj,jl) = 0._wp 1093 END IF 1094 END DO 1095 DO jl = m_index+1, jpl 1096 h_ip(ji,jj,jl) = 0._wp 1097 a_ip(ji,jj,jl) = 0._wp 1098 v_ip(ji,jj,jl) = 0._wp 1099 END DO 1100 1101 ENDIF 1102 1103 END_2D 1104 1105 END SUBROUTINE ice_thd_pnd_area 1106 1107 1108 SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 1109 !!------------------------------------------------------------------- 1110 !! *** ROUTINE ice_thd_pnd_depth *** 1111 !! 1112 !! ** Purpose : Compute melt pond depth 1113 !!------------------------------------------------------------------- 1114 1115 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1116 aicen, & 1117 asnon, & 1118 hsnon, & 1119 alfan, & 1120 cum_max_vol 1121 1122 REAL (wp), INTENT(IN) :: & 1123 zvolp 1124 1125 REAL (wp), INTENT(OUT) :: & 1126 hpond 1127 1128 INTEGER, INTENT(OUT) :: & 1129 m_index 1130 1131 INTEGER :: n, ns 1132 1133 REAL (wp), DIMENSION(0:jpl+1) :: & 1134 hitl, & 1135 aicetl 1136 1137 REAL (wp) :: & 1138 rem_vol, & 1139 area, & 1140 vol, & 1141 tmp, & 1142 z0 = 0.0_wp 1143 1144 !---------------------------------------------------------------- 1145 ! hpond is zero if zvolp is zero - have we fully drained? 1146 !---------------------------------------------------------------- 1147 1148 IF (zvolp < epsi10) THEN 1149 hpond = z0 1150 m_index = 0 1151 ELSE 1152 1153 !---------------------------------------------------------------- 1154 ! Calculate the category where water fills up to 1155 !---------------------------------------------------------------- 1156 1157 !----------| 1158 ! | 1159 ! | 1160 ! |----------| -- -- 1161 !__________|__________|_________________________________________ ^ 1162 ! | | rem_vol ^ | Semi-filled 1163 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1164 ! | | | | 1165 ! | | | |hpond 1166 ! | | |----------| | |------- 1167 ! | | | | | | 1168 ! | | | |---v-----| 1169 ! | | m_index | | | 1170 !------------------------------------------------------------- 1171 1172 m_index = 0 ! 1:m_index categories have water in them 1173 DO n = 1, jpl 1174 IF (zvolp <= cum_max_vol(n)) THEN 1175 m_index = n 1176 IF (n == 1) THEN 1177 rem_vol = zvolp 1178 ELSE 1179 rem_vol = zvolp - cum_max_vol(n-1) 1180 END IF 1181 exit ! to break out of the loop 1182 END IF 1183 END DO 1184 m_index = min(jpl-1, m_index) 1185 1186 !---------------------------------------------------------------- 1187 ! semi-filled layer may have m_index different snow in it 1188 !---------------------------------------------------------------- 1189 1190 !----------------------------------------------------------- ^ 1191 ! | alfan(m_index+1) 1192 ! | 1193 !hitl(3)--> |----------| | 1194 !hitl(2)--> |------------| * * * * *| | 1195 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1196 !hitl(0)-->------------------------------------------------- | ^ 1197 ! various snow from lower categories | |alfa(m_index) 1198 1199 ! hitl - heights of the snow layers from thinner and current categories 1200 ! aicetl - area of each snow depth in this layer 1201 1202 hitl(:) = z0 1203 aicetl(:) = z0 1204 DO n = 1, m_index 1205 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1206 alfan(m_index+1) - alfan(m_index)), z0) 1207 aicetl(n) = asnon(n) 1208 1209 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1210 END DO 1211 1212 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1213 aicetl(m_index+1) = z0 1214 1215 !---------------------------------------------------------------- 1216 ! reorder array according to hitl 1217 ! snow heights not necessarily in height order 1218 !---------------------------------------------------------------- 1219 1220 DO ns = 1, m_index+1 1221 DO n = 0, m_index - ns + 1 1222 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1223 tmp = hitl(n) 1224 hitl(n) = hitl(n+1) 1225 hitl(n+1) = tmp 1226 tmp = aicetl(n) 1227 aicetl(n) = aicetl(n+1) 1228 aicetl(n+1) = tmp 1229 END IF 1230 END DO 1231 END DO 1232 1233 !---------------------------------------------------------------- 1234 ! divide semi-filled layer into set of sublayers each vertically homogenous 1235 !---------------------------------------------------------------- 1236 1237 !hitl(3)---------------------------------------------------------------- 1238 ! | * * * * * * * * 1239 ! |* * * * * * * * * 1240 !hitl(2)---------------------------------------------------------------- 1241 ! | * * * * * * * * | * * * * * * * * 1242 ! |* * * * * * * * * |* * * * * * * * * 1243 !hitl(1)---------------------------------------------------------------- 1244 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1245 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1246 !hitl(0)---------------------------------------------------------------- 1247 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1248 1249 ! move up over layers incrementing volume 1250 DO n = 1, m_index+1 1251 1252 area = sum(aicetl(:)) - & ! total area of sub-layer 1253 (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1254 1255 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1256 1257 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1258 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1259 1260 exit 1261 ELSE ! still in sub-layer below the sub-layer with the depth 1262 rem_vol = rem_vol - vol 1263 END IF 1264 1265 END DO 1266 1267 END IF 1268 1269 END SUBROUTINE ice_thd_pnd_depth 1270 1271 1272 SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 1273 !!------------------------------------------------------------------- 1274 !! *** ROUTINE ice_thd_pnd_perm *** 1275 !! 1276 !! ** Purpose : Determine the liquid fraction of brine in the ice 1277 !! and its permeability 1278 !!------------------------------------------------------------------- 1279 1280 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1281 ticen, & ! internal ice temperature (K) 1282 salin ! salinity (ppt) !js: ppt according to cice 1283 1284 REAL (wp), INTENT(OUT) :: & 1285 perm ! permeability 1286 1287 REAL (wp) :: & 1288 Sbr ! brine salinity 1289 1290 REAL (wp), DIMENSION(nlay_i) :: & 1291 Tin, & ! ice temperature 1292 phi ! liquid fraction 1293 1294 INTEGER :: k 1295 1296 !----------------------------------------------------------------- 1297 ! Compute ice temperatures from enthalpies using quadratic formula 1298 !----------------------------------------------------------------- 1299 1300 DO k = 1,nlay_i 1301 Tin(k) = ticen(k) - rt0 !js: from K to degC 1302 END DO 1303 1304 !----------------------------------------------------------------- 1305 ! brine salinity and liquid fraction 1306 !----------------------------------------------------------------- 1307 1308 DO k = 1, nlay_i 1309 1310 Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1311 ! Best expression to date is that one (Vancoppenolle et al JGR 2019) 1312 ! Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1313 phi(k) = salin(k) / Sbr 1314 1315 END DO 1316 1317 !----------------------------------------------------------------- 1318 ! permeability 1319 !----------------------------------------------------------------- 1320 1321 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1322 1323 END SUBROUTINE ice_thd_pnd_perm 1324 477 1325 SUBROUTINE ice_thd_pnd_init 478 1326 !!------------------------------------------------------------------- … … 511 1359 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 512 1360 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 513 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flus 1361 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush 514 1362 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 515 1363 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
Note: See TracChangeset
for help on using the changeset viewer.
