Changeset 14669
- Timestamp:
- 2021-04-01T11:55:21+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_ponds/src/ICE/icethd_pnd.F90
r14659 r14669 520 520 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 521 521 zTp, & ! pond freezing temperature (C) 522 lid_Ttop,& ! lid top temperature (K) 522 523 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 523 524 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 524 525 z1_rhow, & ! inverse water density 525 zv_pnd , & ! volume of meltwater contributing to ponds526 526 zv_mlt ! total amount of meltwater produced 527 527 … … 530 530 zvolp_res !! remaining melt pond water available after drainage 531 531 532 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 532 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, & 533 zv_pnd !! volume of meltwater contributing to ponds 533 534 534 535 INTEGER :: ji, jj, jk, jl ! loop indices … … 594 595 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 595 596 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 596 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl597 zv_pnd(ji,jj,jl) = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 597 598 598 599 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice ! diags … … 608 609 609 610 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 610 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness611 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd(ji,jj,jl) ! use pond water to increase thickness 611 612 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 612 613 … … 652 653 DO jl = 1, jpl-1 653 654 654 IF ( v_il(ji,jj,jl) > epsi10 ) THEN655 655 656 !---------------------------------------------------------------- 657 ! Lid melting: floating upper ice layer melts in whole or part 658 !---------------------------------------------------------------- 659 ! Use Tsfc for each category 660 IF ( t_su(ji,jj,jl) > zTp ) THEN 661 662 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 663 664 IF ( zdvice > epsi10 ) THEN 665 666 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 667 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 668 ! as it is already counted in surface melt 669 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 670 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 671 672 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 673 ! ice lid melted and category is pond covered 674 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 675 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 676 v_il(ji,jj,jl) = 0._wp 677 ENDIF 678 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 679 680 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 681 682 ENDIF 683 684 !---------------------------------------------------------------- 685 ! Freeze pre-existing lid 686 !---------------------------------------------------------------- 687 688 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 689 690 ! differential growth of base of surface floating ice layer 691 ! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 ! line valid for temp in C 692 zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0 693 zomega = rcnd_i * zdTice / zrhoi_L 694 zdHui = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 695 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 696 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 656 !---------------------------------------------------------------- 657 ! Lid melting: floating upper ice layer melts in whole or part 658 !---------------------------------------------------------------- 659 660 zdvice = MIN( zv_pnd(ji,jj,jl), v_il(ji,jj,jl) ) 661 662 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 663 664 IF ( v_il(ji,jj,jl) < epsi10 ) THEN 665 ! ice lid melted 666 v_il(ji,jj,jl) = 0._wp 667 ENDIF 668 669 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 670 671 !---------------------------------------------------------------- 672 ! Freeze lid 673 !---------------------------------------------------------------- 674 675 ! The temperature of the lid's top surface is warmer than the ice top temperature 676 ! by about 2 degrees (warmed by the melt pond under the lid) 677 lid_Ttop = t_su(ji,jj,jl) + 2.0_wp 678 679 ! differential growth of base of surface floating ice layer 680 zdTice = MAX( - ( lid_Ttop - zTp ) , 0._wp ) ! > 0 681 zomega = rcnd_i * zdTice / zrhoi_L 682 zdvice = SQRT( 2._wp * zomega * rdt_ice * a_ip(ji,jj,jl)**2 + v_il(ji,jj,jl)**2 ) & 683 - v_il(ji,jj,jl) 684 zdvice = min( zdvice , v_ip(ji,jj,jl) ) 685 686 IF ( zdvice > epsi10 ) THEN 687 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 688 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 689 690 ! Possibly this check might not be needed to avoid divide by zero errors as at 691 ! the end of ice_thd_pnd_area all ponds with an area less than epsi10 were zerod anyway 692 IF ( a_ip(ji,jj,jl) > epsi10 ) THEN 693 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 694 ENDIF 695 696 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 697 698 ENDIF 697 699 698 IF ( zdvice > epsi10 ) THEN699 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice700 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice701 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice702 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice703 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)704 705 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag706 707 ENDIF708 709 ENDIF ! Tsfcn(i,j,n)710 711 700 !---------------------------------------------------------------- 712 ! Freeze new lids701 ! Remove ponds if lids are much larger than ponds 713 702 !---------------------------------------------------------------- 714 ! upper ice layer begins to form 715 ! note: albedo does not change 716 717 ELSE ! v_il < epsi10 718 719 ! thickness of newly formed ice 720 ! the surface temperature of a meltpond is the same as that 721 ! of the ice underneath (0C), and the thermodynamic surface 722 ! flux is the same 723 724 !!! we need net surface energy flux, excluding conduction 725 !!! fsurf is summed over categories in CICE 726 !!! we have the category-dependent flux, let us use it ? 727 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 728 zdHui = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp ) 729 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 730 731 IF ( zdvice > epsi10 ) THEN 732 733 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 734 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 735 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 736 737 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 738 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 739 740 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 703 704 IF ( a_ip(ji,jj,jl) > epsi10) THEN 705 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) 706 IF ( h_il(ji,jj,jl) > h_ip(ji,jj,jl) * 10._wp ) THEN 707 a_ip(ji,jj,jl) = 0._wp 708 h_ip(ji,jj,jl) = 0._wp 709 v_ip(ji,jj,jl) = 0._wp 710 h_il(ji,jj,jl) = 0._wp 711 v_il(ji,jj,jl) = 0._wp 741 712 ENDIF 742 743 ENDIF ! v_il 713 ENDIF 744 714 745 715 END DO ! jl
Note: See TracChangeset
for help on using the changeset viewer.
