Changeset 13879
- Timestamp:
- 2020-11-26T08:31:17+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
r13850 r13879 40 40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 41 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 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) 42 47 43 48 !! * Substitutions … … 441 446 zTp, & ! pond freezing temperature (C) 442 447 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 443 zdvn , & ! change in melt pond volume for fresh water budget (not used for now)444 448 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 445 449 z1_rhow, & ! inverse water density … … 447 451 zdum ! dummy variable 448 452 449 REAL(wp), PARAMETER :: & 450 zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m) 451 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 452 zvp_min = 1.e-4_wp ! minimum pond volume (m) 453 454 REAL(wp), DIMENSION(jpi,jpj) :: zvolp !! total melt pond water available before redistribution and drainage 453 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 454 zvolp_res 455 455 456 456 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i … … 466 466 467 467 !--------------------------------------------------------------- 468 ! Initiali ze468 ! Initialise 469 469 !--------------------------------------------------------------- 470 470 … … 487 487 END WHERE 488 488 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 489 490 !----------------------------------------------------------------- 491 ! Start pond loop 492 !----------------------------------------------------------------- 489 490 !--------------------------------------------------------------- 491 ! Change 2D to 1D 492 !--------------------------------------------------------------- 493 494 !-------------------------------------------------------------- 495 ! Collect total available pond water 496 !-------------------------------------------------------------- 493 497 494 498 zvolp(:,:) = 0._wp 495 499 496 DO jj = 1, jpj 497 DO ji = 1, jpi 498 499 !-------------------------------------------------------------- 500 ! Collect meltwater on ponds 501 !-------------------------------------------------------------- 502 DO jl = 1, jpl 503 500 DO jl = 1, jpl 501 DO jj = 1, jpj 502 DO ji = 1, jpi 503 504 504 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 505 505 … … 518 518 !--- Deepen existing ponds before redistribution and drainage(later on) 519 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) 520 zvolp(ji,jj) = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 521 521 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless 522 522 523 523 ENDIF ! a_i 524 524 525 END DO 526 527 h_ip(ji,jj,:) = 0._wp ! pond thickness 528 a_ip(ji,jj,:) = 0._wp ! pond fraction per grid area 529 530 !----------------------------------------------------------------- 531 ! Identify grid cells with ice 532 !----------------------------------------------------------------- 525 END DO! jl 526 END DO ! jj 527 END DO ! ji 528 529 h_ip(:,:,:) = 0._wp ! pond thickness 530 a_ip(:,:,:) = 0._wp ! pond fraction per grid area 531 532 !-------------------------------------------------------------- 533 ! Redistribute and drain water from ponds 534 !-------------------------------------------------------------- 535 CALL ice_thd_pnd_area( zvolp, zvolp_res ) 536 537 !-------------------------------------------------------------- 538 ! Freeze and melt lid 539 !-------------------------------------------------------------- 540 DO jj = 1, jpj 541 DO ji = 1, jpi 533 542 534 543 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 535 536 !---------------------------------------------------------------------------- 537 ! Calculate pond area and depth (redistribution of melt water on topography) 538 !---------------------------------------------------------------------------- 539 ! CALL ice_thd_pnd_area( at_i(ji,jj) , vt_i(ji,jj), & 540 ! a_i (ji,jj,:), v_i(ji,jj,:), v_s(ji,jj,:), & 541 ! t_i(ji,jj,:,:), sz_i(ji,jj,:,:), & 542 ! v_ip(ji,jj,:) , zvolp(ji,jj), & 543 ! a_ip(ji,jj,:), h_ip(ji,jj,:), zdvn) 544 545 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvn ! ---> zdvn is drainage (useless for now since we dont change fw budget, no ?) 546 544 547 545 !-------------------------- 548 546 ! Pond lid growth and melt … … 565 563 IF ( t_su(ji,jj,jl) > zTp ) THEN 566 564 567 zdvice = min( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl))565 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 568 566 IF ( zdvice > epsi10 ) then 569 567 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice … … 617 615 ! flux is the same 618 616 619 !!!!! FSURF !!!!!620 617 !!! we need net surface energy flux, excluding conduction 621 618 !!! fsurf is summed over categories in CICE … … 639 636 640 637 v_ip(ji,jj,:) = 0._wp 641 v_il (ji,jj,:)= 0._wp638 v_il(ji,jj,:) = 0._wp 642 639 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 643 640 ! zvolp(ji,jj) = 0._wp … … 645 642 ENDIF 646 643 644 END DO ! jj 645 END DO ! ji 646 647 647 !--------------------------------------------------------------- 648 ! remove ice lid if there is no liquid pond 649 ! v_il may be nonzero on category ncat due to dynamics 648 ! Clean-up variables (probably duplicates what icevar would do) 650 649 !--------------------------------------------------------------- 650 ! MV comment 651 ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 652 ! icevar should recompute all other variables (if needed at all) 651 653 652 654 DO jl = 1, jpl 655 DO jj = 1, jpj 656 DO ji = 1, jpi 657 658 IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 659 .AND. v_il (ji,jj,jl) > epsi10) THEN 660 v_il(ji,jj,jl) = 0._wp 661 ENDIF 653 662 654 IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 655 .AND. v_il (ji,jj,jl) > epsi10) THEN 656 v_il(ji,jj,jl) = 0._wp 657 ENDIF 658 659 ! reload tracers 660 IF ( a_ip(ji,jj,jl) > epsi10) THEN 661 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 662 ELSE 663 v_il(ji,jj,jl) = 0._wp 664 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar 665 ENDIF 666 IF ( a_ip(ji,jj,jl) > epsi10 ) THEN 667 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 668 !h_ip(ji,jj,jl) = zhpondn(ji,jj,jl) 669 ELSE 670 a_ip_frac(ji,jj,jl) = 0._wp 671 h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 672 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 673 ENDIF 663 ! reload tracers 664 IF ( a_ip(ji,jj,jl) > epsi10) THEN 665 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 666 ELSE 667 v_il(ji,jj,jl) = 0._wp 668 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar 669 ENDIF 670 671 IF ( a_ip(ji,jj,jl) > epsi10 ) THEN 672 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 673 !h_ip(ji,jj,jl) = zhpondn(ji,jj,jl) 674 ELSE 675 a_ip_frac(ji,jj,jl) = 0._wp 676 h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 677 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 678 ENDIF 674 679 675 END DO ! jl 676 677 END DO ! jj 678 END DO ! ji 680 END DO ! ji 681 END DO ! jj 682 END DO ! jl 679 683 680 684 END SUBROUTINE pnd_TOPO 681 685 682 686 683 SUBROUTINE ice_thd_pnd_area(aice, vice, & 684 aicen, vicen, vsnon, ticen, & 685 salin, zvolpn, zvolp, & 686 zapondn,zhpondn,dvolp) 687 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 687 688 688 689 !!------------------------------------------------------------------- 689 690 !! *** ROUTINE ice_thd_pnd_area *** 690 691 !! 691 !! ** Purpose : Given the total volume of meltwater, update692 !! pond fraction (a_ip) and depth (should be volume)692 !! ** Purpose : Given the total volume of available pond water, 693 !! redistribute and drain water 693 694 !! 694 695 !! ** 695 696 !! 696 697 !!------------------------------------------------------------------ 697 698 REAL (wp), INTENT(IN) :: & 699 aice, vice 700 701 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 702 aicen, vicen, vsnon 703 704 REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 705 ticen, salin 706 707 REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 708 zvolpn 709 710 REAL (wp), INTENT(INOUT) :: & 711 zvolp, dvolp 712 713 REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 714 zapondn, zhpondn 698 699 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 700 zvolp, & ! total available pond water 701 zdvolp ! remaining meltwater after redistribution 715 702 716 703 INTEGER :: & 717 n , ns, &704 ns, & 718 705 m_index, & 719 706 permflag … … 744 731 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 745 732 733 INTEGER :: ji, jj, jk, jl ! loop indices 734 746 735 !-----------| 747 736 ! | … … 752 741 ! | | | | 753 742 ! | | | |-----------| |------- 754 ! | | |alfan( n)| | |743 ! | | |alfan(jl)| | | 755 744 ! | | | | |--------------| 756 745 ! | | | | | | … … 758 747 ! | | ^ | | | 759 748 ! | | | | |--------------| 760 ! | | |betan( n)| | |749 ! | | |betan(jl)| | | 761 750 ! | | | |-----------| |------- 762 751 ! | | | | … … 767 756 !-----------| 768 757 758 a_ip(:,:,:) = 0._wp 759 h_ip(:,:,:) = 0._wp 760 761 DO jj = 1, jpj 762 DO ji = 1, jpi 763 764 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 765 769 766 !------------------------------------------------------------------- 770 767 ! initialize 771 768 !------------------------------------------------------------------- 772 769 773 DO n= 1, jpl774 775 zapondn(n) = 0._wp776 zhpondn(n) = 0._wp770 DO jl = 1, jpl 771 772 a_ip(ji,jj,jl) = 0._wp 773 h_ip(ji,jj,jl) = 0._wp 777 774 778 775 !---------------------------------------- 779 ! X)compute the effective snow fraction776 ! compute the effective snow fraction 780 777 !---------------------------------------- 781 IF (aicen(n) < epsi10) THEN 782 hicen(n) = 0._wp 783 hsnon(n) = 0._wp 784 reduced_aicen(n) = 0._wp 785 asnon(n) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 778 779 IF (a_i(ji,jj,jl) < epsi10) THEN 780 hicen(jl) = 0._wp 781 hsnon(jl) = 0._wp 782 reduced_aicen(jl) = 0._wp 783 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 786 784 ELSE 787 hicen( n) = vicen(n) / aicen(n)788 hsnon( n) = vsnon(n) / aicen(n)789 reduced_aicen( n) = 1._wp ! n=jpl785 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 786 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 787 reduced_aicen(jl) = 1._wp ! n=jpl 790 788 791 789 !js: initial code in NEMO_DEV 792 !IF (n < jpl) reduced_aicen( n) = aicen(n) &793 ! * (-0.024_wp*hicen( n) + 0.832_wp)790 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 791 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 794 792 795 793 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 796 IF ( n < jpl) reduced_aicen(n) = aicen(n) &797 * max(0.2_wp,(-0.024_wp*hicen( n) + 0.832_wp))798 799 asnon( n) = reduced_aicen(n) ! effective snow fraction (empirical)794 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 795 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 796 797 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 800 798 ! MV should check whether this makes sense to have the same effective snow fraction in here 801 799 ! OLI: it probably doesn't … … 817 815 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 818 816 819 alfan( n) = 0.6 * hicen(n)820 betan( n) = 0.4 * hicen(n)821 822 cum_max_vol( n) = 0._wp823 cum_max_vol_tmp( n) = 0._wp817 alfan(jl) = 0.6 * hicen(jl) 818 betan(jl) = 0.4 * hicen(jl) 819 820 cum_max_vol(jl) = 0._wp 821 cum_max_vol_tmp(jl) = 0._wp 824 822 825 823 END DO ! jpl … … 827 825 cum_max_vol_tmp(0) = 0._wp 828 826 drain = 0._wp 829 dvolp= 0._wp827 zdvolp(ji,jj) = 0._wp 830 828 831 829 !---------------------------------------------------------- 832 ! x)Drain overflow water, update pond fraction and volume830 ! Drain overflow water, update pond fraction and volume 833 831 !---------------------------------------------------------- 834 832 … … 836 834 ! the maximum amount of water that can be contained up to each ice category 837 835 !-------------------------------------------------------------------------- 838 839 ! MV840 836 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 841 837 ! Then the excess volume cum_max_vol(jl) drains out of the system 842 838 ! It should be added to wfx_pnd_out 843 ! END MV 844 !js 18/04/19: XXX do something about this flux thing 845 846 DO n = 1, jpl-1 ! last category can not hold any volume 847 848 IF (alfan(n+1) >= alfan(n) .AND. alfan(n+1) > 0._wp ) THEN 839 840 DO jl = 1, jpl-1 ! last category can not hold any volume 841 842 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 849 843 850 844 ! total volume in level including snow 851 cum_max_vol_tmp( n) = cum_max_vol_tmp(n-1) + &852 (alfan( n+1) - alfan(n)) * sum(reduced_aicen(1:n))845 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 846 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 853 847 854 848 ! subtract snow solid volumes from lower categories in current level 855 DO ns = 1, n856 cum_max_vol_tmp( n) = cum_max_vol_tmp(n) &849 DO ns = 1, jl 850 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 857 851 - rhos/rhow * & ! free air fraction that can be filled by water 858 852 asnon(ns) * & ! effective areal fraction of snow in that category 859 max(min(hsnon(ns)+alfan(ns)-alfan( n), alfan(n+1)-alfan(n)), 0._wp)853 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 860 854 END DO 861 855 862 856 ELSE ! assume higher categories unoccupied 863 cum_max_vol_tmp( n) = cum_max_vol_tmp(n-1)857 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 864 858 END IF 865 !IF (cum_max_vol_tmp( n) < z0) THEN859 !IF (cum_max_vol_tmp(jl) < z0) THEN 866 860 ! CALL abort_ice('negative melt pond volume') 867 861 !END IF … … 873 867 ! is there more meltwater than can be held in the floe? 874 868 !---------------------------------------------------------------- 875 IF (zvolp >= cum_max_vol(jpl)) THEN876 drain = zvolp - cum_max_vol(jpl) + epsi10877 zvolp = zvolp- drain ! update meltwater volume available878 dvolp= drain ! this is the drained water879 IF (zvolp < epsi10) THEN880 dvolp = dvolp + zvolp881 zvolp = 0._wp869 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 870 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 871 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 872 zdvolp(ji,jj) = drain ! this is the drained water 873 IF (zvolp(ji,jj) < epsi10) THEN 874 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 875 zvolp(ji,jj) = 0._wp 882 876 END IF 883 877 END IF 884 878 885 879 ! height and area corresponding to the remaining volume 886 887 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 888 889 DO n=1, m_index 890 !zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 880 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 881 882 DO jl = 1, m_index 883 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 891 884 ! ! volume instead, no ? 892 zhpondn(n) = max((hpond - alfan(n) + alfan(1)), 0._wp) !js: from CICE 5.1.2893 zapondn(n) = reduced_aicen(n)885 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 886 a_ip(ji,jj,jl) = reduced_aicen(jl) 894 887 ! in practise, pond fraction depends on the empirical snow fraction 895 888 ! so in turn on ice thickness 896 889 END DO 897 !zapond = sum( zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d890 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 898 891 899 892 !------------------------------------------------------------------------ … … 903 896 904 897 ! sea water level 905 msno = 0._wp906 DO n=1,jpl907 msno = msno + v snon(n) * rhos898 msno = 0._wp 899 DO jl = 1 , jpl 900 msno = msno + v_s(ji,jj,jl) * rhos 908 901 END DO 909 floe_weight = ( msno + rhoi*vice + rau0*zvolp) / aice902 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 910 903 hsl_rel = floe_weight / rau0 & 911 - ( (sum(betan(:)*aicen(:))/aice) + alfan(1))904 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 912 905 913 906 deltah = hpond - hsl_rel … … 916 909 ! drain if ice is permeable 917 910 permflag = 0 911 918 912 IF (pressure_head > 0._wp) THEN 919 DO n = 1, jpl-1 920 IF (hicen(n) /= 0._wp) THEN 921 !IF (hicen(n) > 0._wp) THEN !js: from CICE 5.1.2 913 DO jl = 1, jpl-1 914 IF ( hicen(jl) /= 0._wp ) THEN 915 916 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 917 922 918 perm = 0._wp ! MV ugly dummy patch 923 CALL ice_thd_pnd_perm(ticen(:,n), salin(:,n), perm)919 ! CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 924 920 IF (perm > 0._wp) permflag = 1 925 921 926 drain = perm* zapondn(n)*pressure_head*rdt_ice / &927 (viscosity*hicen( n))928 dvolp = dvolp + min(drain, zvolp)929 zvolp = max(zvolp- drain, 0._wp)930 IF (zvolp < epsi10) THEN931 dvolp = dvolp + zvolp932 zvolp = 0._wp922 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 923 (viscosity*hicen(jl)) 924 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 925 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 926 IF (zvolp(ji,jj) < epsi10) THEN 927 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 928 zvolp(ji,jj) = 0._wp 933 929 END IF 934 930 END IF … … 938 934 IF (permflag > 0) THEN 939 935 ! recompute pond depth 940 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp , cum_max_vol, hpond, m_index)941 DO n=1, m_index942 zhpondn(n) = hpond - alfan(n) + alfan(1)943 zapondn(n) = reduced_aicen(n)936 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 937 DO jl = 1, m_index 938 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 939 a_ip(ji,jj,jl) = reduced_aicen(jl) 944 940 END DO 945 !zapond = sum( zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d941 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 946 942 END IF 947 943 END IF ! pressure_head 948 944 949 945 !------------------------------- 950 ! X)remove water from the snow946 ! remove water from the snow 951 947 !------------------------------- 952 948 !------------------------------------------------------------------------ … … 956 952 957 953 ! Calculate pond volume for lower categories 958 DO n=1,m_index-1959 zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow960 - (rhos/rhow) * asnon(n) * min(hsnon(n), zhpondn(n))954 DO jl = 1,m_index-1 955 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)) 961 957 END DO 962 958 … … 968 964 ! m_index = last category with melt pond 969 965 970 IF (m_index == 1) zvolpn(m_index) = zvolp! volume of mw in 1st category is the total volume of melt water966 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 971 967 972 968 IF (m_index > 1) THEN 973 IF (zvolp > sum(zvolpn(1:m_index-1))) THEN974 zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1))969 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 970 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 975 971 ELSE 976 zvolpn(m_index) = 0._wp977 zhpondn(m_index) = 0._wp978 zapondn(m_index) = 0._wp979 ! If remaining pond volume is negative reduce pond volume of980 ! lower category981 IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &982 zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp972 v_ip(ji,jj,m_index) = 0._wp 973 h_ip(ji,jj,m_index) = 0._wp 974 a_ip(ji,jj,m_index) = 0._wp 975 ! If remaining pond volume is negative reduce pond volume of 976 ! lower category 977 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 978 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) 983 979 END IF 984 980 END IF 985 981 986 DO n=1,m_index987 IF ( zapondn(n) > epsi10) THEN988 zhpondn(n) = zvolpn(n) / zapondn(n)982 DO jl = 1,m_index 983 IF (a_ip(ji,jj,jl) > epsi10) THEN 984 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 989 985 ELSE 990 dvolp = dvolp + zvolpn(n)991 zhpondn(n) = 0._wp992 zvolpn(n) = 0._wp993 zapondn(n) = 0._wp986 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 987 h_ip(ji,jj,jl) = 0._wp 988 v_ip(ji,jj,jl) = 0._wp 989 a_ip(ji,jj,jl) = 0._wp 994 990 END IF 995 991 END DO 996 DO n= m_index+1, jpl997 zhpondn(n) = 0._wp998 zapondn(n) = 0._wp999 zvolpn (n) = 0._wp992 DO jl = m_index+1, jpl 993 h_ip(ji,jj,jl) = 0._wp 994 a_ip(ji,jj,jl) = 0._wp 995 v_ip(ji,jj,jl) = 0._wp 1000 996 END DO 997 998 ENDIF 999 END DO ! ji 1000 END DO ! jj 1001 1001 1002 1002 END SUBROUTINE ice_thd_pnd_area
Note: See TracChangeset
for help on using the changeset viewer.
