Changeset 14072 for NEMO/trunk/src/ICE/icethd_pnd.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icethd_pnd.F90 (modified) (67 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icethd_pnd.F90
r14005 r14072 1 MODULE icethd_pnd 1 MODULE icethd_pnd 2 2 !!====================================================================== 3 3 !! *** MODULE icethd_pnd *** … … 41 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 42 42 43 !-------------------------------------------------------------------------- 43 !-------------------------------------------------------------------------- 44 44 ! Diagnostics for pond volume per area 45 45 ! … … 56 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_drn ! pond volume lost by drainage [-] 57 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_lid ! exchange with lid / refreezing [-] 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] 59 59 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_mlt_1d ! meltwater pond volume input [-] 60 60 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_drn_1d ! pond volume lost by drainage [-] … … 75 75 !!------------------------------------------------------------------- 76 76 !! *** ROUTINE ice_thd_pnd *** 77 !! 77 !! 78 78 !! ** Purpose : change melt pond fraction and thickness 79 79 !! … … 84 84 INTEGER :: ji, jj, jl ! loop indices 85 85 !!------------------------------------------------------------------- 86 86 87 87 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 88 88 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) … … 111 111 END_2D 112 112 END DO 113 113 114 114 !------------------------------ 115 115 ! Identify grid cells with ice … … 148 148 DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf ) 149 149 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 150 151 END SUBROUTINE ice_thd_pnd 152 153 154 SUBROUTINE pnd_CST 150 151 END SUBROUTINE ice_thd_pnd 152 153 154 SUBROUTINE pnd_CST 155 155 !!------------------------------------------------------------------- 156 156 !! *** ROUTINE pnd_CST *** … … 158 158 !! ** Purpose : Compute melt pond evolution 159 159 !! 160 !! ** Method : Melt pond fraction and thickness are prescribed 160 !! ** Method : Melt pond fraction and thickness are prescribed 161 161 !! to non-zero values when t_su = 0C 162 162 !! 163 163 !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) 164 !! 164 !! 165 165 !! ** Note : Coupling with such melt ponds is only radiative 166 166 !! Advection, ridging, rafting... are bypassed … … 172 172 !!------------------------------------------------------------------- 173 173 DO jl = 1, jpl 174 174 175 175 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 176 176 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) … … 185 185 ! 186 186 IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 187 h_ip_1d(ji) = rn_hpnd 187 h_ip_1d(ji) = rn_hpnd 188 188 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 189 189 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 190 190 ELSE 191 h_ip_1d(ji) = 0._wp 191 h_ip_1d(ji) = 0._wp 192 192 a_ip_1d(ji) = 0._wp 193 193 h_il_1d(ji) = 0._wp … … 222 222 !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing 223 223 !! We work with volumes and then redistribute changes into thickness and concentration 224 !! assuming linear relationship between the two. 224 !! assuming linear relationship between the two. 225 225 !! 226 226 !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- … … 237 237 !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- 238 238 !! 239 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 239 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 240 240 !! H = lid thickness 241 241 !! Lf = latent heat of fusion … … 260 260 !! 261 261 !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 262 !! 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 262 !! 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 264 264 !! 265 265 !! ** References : Flocco and Feltham (JGR, 2007) … … 267 267 !! Holland et al (J. Clim, 2012) 268 268 !! Hunke et al (OM 2012) 269 !!------------------------------------------------------------------- 269 !!------------------------------------------------------------------- 270 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 271 271 !! … … 287 287 INTEGER :: ji, jk, jl ! loop indices 288 288 !!------------------------------------------------------------------- 289 z1_rhow = 1._wp / rhow 289 z1_rhow = 1._wp / rhow 290 290 z1_aspect = 1._wp / zaspect 291 z1_Tp = 1._wp / zTp 292 291 z1_Tp = 1._wp / zTp 292 293 293 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti), at_i ) 294 294 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 295 295 296 296 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 297 297 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) … … 315 315 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 316 316 END DO 317 317 318 318 !----------------------- 319 319 ! Melt pond calculations … … 342 342 zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 343 343 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 344 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 344 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 345 345 ! 346 346 !--- overflow ---! … … 349 349 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 350 350 ! a_ip_max = zfr_mlt * a_i 351 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 351 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 352 352 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 353 353 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) … … 356 356 ! If pond depth exceeds half the ice thickness then reduce the pond volume 357 357 ! h_ip_max = 0.5 * h_i 358 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 358 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 359 359 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 360 360 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) … … 375 375 IF( ln_pnd_lids ) THEN 376 376 ! 377 !--- Lid growing and subsequent pond shrinking ---! 377 !--- Lid growing and subsequent pond shrinking ---! 378 378 zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 379 379 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors … … 386 386 387 387 ELSE 388 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 388 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 389 389 ! Pond shrinking 390 390 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) … … 398 398 ! 399 399 400 !------------------------------------------------! 400 !------------------------------------------------! 401 401 ! Pond drainage through brine network (flushing) ! 402 402 !------------------------------------------------! … … 420 420 ! Do the drainage using Darcy's law 421 421 zdv_flush = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 423 423 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 424 424 … … 479 479 480 480 481 SUBROUTINE pnd_TOPO 482 481 SUBROUTINE pnd_TOPO 482 483 483 !!------------------------------------------------------------------- 484 484 !! *** ROUTINE pnd_TOPO *** … … 488 488 !! 489 489 !! ** Method : This code is initially based on Flocco and Feltham 490 !! (2007) and Flocco et al. (2010). 490 !! (2007) and Flocco et al. (2010). 491 491 !! 492 492 !! - Calculate available pond water base on surface meltwater … … 532 532 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 533 533 zvolp_res !! remaining melt pond water available after drainage 534 534 535 535 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 536 536 … … 545 545 546 546 CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 547 547 548 548 !--------------------------------------------------------------- 549 549 ! Initialise … … 553 553 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 554 554 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 555 z1_rhow = 1._wp / rhow 555 z1_rhow = 1._wp / rhow 556 556 557 557 ! Set required ice variables (hard-coded here for now) 558 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 559 558 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 559 560 560 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 561 561 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 562 562 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 563 563 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 564 564 565 565 ! thickness 566 566 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) … … 568 568 END WHERE 569 569 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 570 570 571 571 !--------------------------------------------------------------- 572 572 ! Change 2D to 1D 573 573 !--------------------------------------------------------------- 574 ! MV 574 ! MV 575 575 ! a less computing-intensive version would have 2D-1D passage here 576 576 ! use what we have in iceitd.F90 (incremental remapping) … … 582 582 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 583 583 ! I cite her words, they are very talkative 584 ! "grid cells with very little ice cover (and hence more open water area) 584 ! "grid cells with very little ice cover (and hence more open water area) 585 585 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 586 586 ! "This results in the same runoff fraction r for each ice category within a grid cell" 587 587 588 588 zvolp(:,:) = 0._wp 589 589 590 590 DO jl = 1, jpl 591 591 DO_2D( 1, 1, 1, 1 ) 592 592 593 593 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 594 594 595 595 !--- Available and contributing meltwater for melt ponding ---! 596 596 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 … … 601 601 602 602 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice ! diags 603 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice 603 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice 604 604 605 605 !--- Create possible new ponds … … 610 610 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 611 611 ENDIF 612 612 613 613 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 614 614 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 615 615 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 616 616 617 617 !--- Total available pond water volume (pre-existing + newly produced)j 618 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 618 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 619 619 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 620 620 621 621 ENDIF ! a_i 622 622 623 623 END_2D 624 624 END DO ! ji 625 625 626 626 !-------------------------------------------------------------- 627 627 ! Redistribute and drain water from ponds 628 !-------------------------------------------------------------- 628 !-------------------------------------------------------------- 629 629 CALL ice_thd_pnd_area( zvolp, zvolp_res ) 630 630 631 631 !-------------------------------------------------------------- 632 632 ! Melt pond lid growth and melt 633 !-------------------------------------------------------------- 634 633 !-------------------------------------------------------------- 634 635 635 IF( ln_pnd_lids ) THEN 636 636 … … 638 638 639 639 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 640 640 641 641 !-------------------------- 642 642 ! Pond lid growth and melt … … 648 648 END DO 649 649 zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 650 650 651 651 DO jl = 1, jpl-1 652 652 653 653 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 654 654 655 655 !---------------------------------------------------------------- 656 656 ! Lid melting: floating upper ice layer melts in whole or part … … 660 660 661 661 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 662 662 663 663 IF ( zdvice > epsi10 ) THEN 664 664 665 665 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 666 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 666 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 667 667 ! as it is already counted in surface melt 668 668 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 669 669 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 670 670 671 671 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 672 672 ! ice lid melted and category is pond covered 673 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 674 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 673 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 674 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 675 675 v_il(ji,jj,jl) = 0._wp 676 676 ENDIF 677 677 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 678 678 679 679 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 680 680 681 681 ENDIF 682 682 683 683 !---------------------------------------------------------------- 684 ! Freeze pre-existing lid 684 ! Freeze pre-existing lid 685 685 !---------------------------------------------------------------- 686 686 … … 688 688 689 689 ! differential growth of base of surface floating ice layer 690 zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 690 zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 691 691 zomega = rcnd_i * zdTice / zrhoi_L 692 692 zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 693 693 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 694 694 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 695 695 696 696 IF ( zdvice > epsi10 ) THEN 697 697 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice … … 700 700 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 701 701 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 702 702 703 703 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 704 704 705 705 ENDIF 706 706 707 707 ENDIF ! Tsfcn(i,j,n) 708 708 … … 714 714 715 715 ELSE ! v_il < epsi10 716 716 717 717 ! thickness of newly formed ice 718 718 ! the surface temperature of a meltpond is the same as that 719 ! of the ice underneath (0C), and the thermodynamic surface 719 ! of the ice underneath (0C), and the thermodynamic surface 720 720 ! flux is the same 721 721 722 722 !!! we need net surface energy flux, excluding conduction 723 723 !!! fsurf is summed over categories in CICE 724 724 !!! we have the category-dependent flux, let us use it ? 725 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 725 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 726 726 zdHui = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) 727 727 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) … … 729 729 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 730 730 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 731 731 732 732 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 733 733 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice … … 735 735 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 736 736 ENDIF 737 737 738 738 ENDIF ! v_il 739 739 740 740 END DO ! jl 741 741 … … 745 745 v_il(ji,jj,:) = 0._wp 746 746 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 747 ! zvolp(ji,jj) = 0._wp 747 ! zvolp(ji,jj) = 0._wp 748 748 749 749 ENDIF … … 769 769 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 770 770 ! ENDIF 771 771 772 772 ! recalculate equivalent pond variables 773 773 IF ( a_ip(ji,jj,jl) > epsi10) THEN … … 779 779 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 780 780 ! ENDIF 781 781 782 782 END_2D 783 783 … … 787 787 END SUBROUTINE pnd_TOPO 788 788 789 789 790 790 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 791 791 … … 793 793 !! *** ROUTINE ice_thd_pnd_area *** 794 794 !! 795 !! ** Purpose : Given the total volume of available pond water, 795 !! ** Purpose : Given the total volume of available pond water, 796 796 !! redistribute and drain water 797 797 !! … … 823 823 !! 824 824 !!------------------------------------------------------------------ 825 825 826 826 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 827 827 zvolp, & ! total available pond water … … 865 865 a_ip(:,:,:) = 0._wp 866 866 h_ip(:,:,:) = 0._wp 867 867 868 868 DO_2D( 1, 1, 1, 1 ) 869 869 870 870 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 871 871 872 872 !------------------------------------------------------------------- 873 873 ! initialize 874 874 !------------------------------------------------------------------- 875 875 876 876 DO jl = 1, jpl 877 877 878 878 !---------------------------------------- 879 879 ! compute the effective snow fraction 880 880 !---------------------------------------- 881 881 882 882 IF (a_i(ji,jj,jl) < epsi10) THEN 883 883 hicen(jl) = 0._wp … … 889 889 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 890 890 reduced_aicen(jl) = 1._wp ! n=jpl 891 891 892 892 !js: initial code in NEMO_DEV 893 893 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 894 894 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 895 895 896 896 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 897 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 897 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 898 898 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 899 899 900 900 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 901 901 ! MV should check whether this makes sense to have the same effective snow fraction in here 902 902 ! OLI: it probably doesn't 903 903 END IF 904 904 905 905 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 906 906 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming … … 911 911 ! alfan = 60% of the ice volume) in each category lies above the reference line, 912 912 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 913 913 914 914 ! MV: 915 915 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 916 916 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 917 917 918 918 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 919 919 920 920 alfan(jl) = 0.6 * hicen(jl) 921 921 betan(jl) = 0.4 * hicen(jl) 922 922 923 923 cum_max_vol(jl) = 0._wp 924 924 cum_max_vol_tmp(jl) = 0._wp 925 925 926 926 END DO ! jpl 927 927 928 928 cum_max_vol_tmp(0) = 0._wp 929 929 drain = 0._wp 930 930 zdvolp(ji,jj) = 0._wp 931 931 932 932 !---------------------------------------------------------- 933 933 ! Drain overflow water, update pond fraction and volume 934 934 !---------------------------------------------------------- 935 935 936 936 !-------------------------------------------------------------------------- 937 937 ! the maximum amount of water that can be contained up to each ice category … … 940 940 ! Then the excess volume cum_max_vol(jl) drains out of the system 941 941 ! It should be added to wfx_pnd_out 942 942 943 943 DO jl = 1, jpl-1 ! last category can not hold any volume 944 944 945 945 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 946 946 947 947 ! total volume in level including snow 948 948 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 949 949 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 950 950 951 951 ! subtract snow solid volumes from lower categories in current level 952 952 DO ns = 1, jl … … 956 956 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 957 957 END DO 958 958 959 959 ELSE ! assume higher categories unoccupied 960 960 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) … … 966 966 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 967 967 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 968 968 969 969 !---------------------------------------------------------------- 970 970 ! is there more meltwater than can be held in the floe? … … 973 973 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 974 974 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 975 975 976 976 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 977 977 978 978 zdvolp(ji,jj) = drain ! this is the drained water 979 979 IF (zvolp(ji,jj) < epsi10) THEN … … 982 982 END IF 983 983 END IF 984 984 985 985 ! height and area corresponding to the remaining volume 986 986 ! routine leaves zvolp unchanged 987 987 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 988 988 989 989 DO jl = 1, m_index 990 990 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update … … 996 996 END DO 997 997 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 998 998 999 999 !------------------------------------------------------------------------ 1000 1000 ! Drainage through brine network (permeability) 1001 1001 !------------------------------------------------------------------------ 1002 1002 !!! drainage due to ice permeability - Darcy's law 1003 1003 1004 1004 ! sea water level 1005 msno = 0._wp 1005 msno = 0._wp 1006 1006 DO jl = 1 , jpl 1007 1007 msno = msno + v_s(ji,jj,jl) * rhos … … 1010 1010 hsl_rel = floe_weight / rho0 & 1011 1011 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1012 1012 1013 1013 deltah = hpond - hsl_rel 1014 1014 pressure_head = grav * rho0 * max(deltah, 0._wp) 1015 1015 1016 1016 ! drain if ice is permeable 1017 1017 permflag = 0 1018 1018 1019 1019 IF (pressure_head > 0._wp) THEN 1020 1020 DO jl = 1, jpl-1 1021 1021 IF ( hicen(jl) /= 0._wp ) THEN 1022 1022 1023 1023 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1024 1024 1025 1025 perm = 0._wp ! MV ugly dummy patch 1026 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1026 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1027 1027 IF (perm > 0._wp) permflag = 1 1028 1028 1029 1029 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & 1030 1030 (viscosity*hicen(jl)) 1031 1031 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1032 1032 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1033 1033 1034 1034 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1035 1035 1036 1036 IF (zvolp(ji,jj) < epsi10) THEN 1037 1037 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) … … 1040 1040 END IF 1041 1041 END DO 1042 1042 1043 1043 ! adjust melt pond dimensions 1044 1044 IF (permflag > 0) THEN … … 1052 1052 END IF 1053 1053 END IF ! pressure_head 1054 1054 1055 1055 !------------------------------- 1056 1056 ! remove water from the snow … … 1060 1060 ! snow in melt ponds is not melted 1061 1061 !------------------------------------------------------------------------ 1062 1062 1063 1063 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1064 1064 ! how much, so I did not diagnose it 1065 1065 ! so if there is a problem here, nobody is going to see it... 1066 1067 1066 1067 1068 1068 ! Calculate pond volume for lower categories 1069 1069 DO jl = 1,m_index-1 … … 1071 1071 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1072 1072 END DO 1073 1073 1074 1074 ! Calculate pond volume for highest category = remaining pond volume 1075 1075 1076 1076 ! The following is completely unclear to Martin at least 1077 1077 ! Could we redefine properly and recode in a more readable way ? 1078 1078 1079 1079 ! m_index = last category with melt pond 1080 1080 1081 1081 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 1082 1082 1083 1083 IF (m_index > 1) THEN 1084 1084 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1085 1085 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1086 1086 ELSE 1087 v_ip(ji,jj,m_index) = 0._wp 1087 v_ip(ji,jj,m_index) = 0._wp 1088 1088 h_ip(ji,jj,m_index) = 0._wp 1089 1089 a_ip(ji,jj,m_index) = 0._wp … … 1094 1094 END IF 1095 1095 END IF 1096 1096 1097 1097 DO jl = 1,m_index 1098 1098 IF (a_ip(ji,jj,jl) > epsi10) THEN … … 1100 1100 ELSE 1101 1101 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1102 h_ip(ji,jj,jl) = 0._wp 1102 h_ip(ji,jj,jl) = 0._wp 1103 1103 v_ip(ji,jj,jl) = 0._wp 1104 1104 a_ip(ji,jj,jl) = 0._wp … … 1106 1106 END DO 1107 1107 DO jl = m_index+1, jpl 1108 h_ip(ji,jj,jl) = 0._wp 1109 a_ip(ji,jj,jl) = 0._wp 1110 v_ip(ji,jj,jl) = 0._wp 1108 h_ip(ji,jj,jl) = 0._wp 1109 a_ip(ji,jj,jl) = 0._wp 1110 v_ip(ji,jj,jl) = 0._wp 1111 1111 END DO 1112 1112 1113 1113 ENDIF 1114 1114 … … 1319 1319 1320 1320 DO k = 1, nlay_i 1321 1321 1322 1322 Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1323 1323 ! Best expression to date is that one (Vancoppenolle et al JGR 2019) 1324 1324 ! Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1325 1325 phi(k) = salin(k) / Sbr 1326 1326 1327 1327 END DO 1328 1328 … … 1335 1335 END SUBROUTINE ice_thd_pnd_perm 1336 1336 1337 SUBROUTINE ice_thd_pnd_init 1337 SUBROUTINE ice_thd_pnd_init 1338 1338 !!------------------------------------------------------------------- 1339 1339 !! *** ROUTINE ice_thd_pnd_init *** … … 1342 1342 !! over sea ice 1343 1343 !! 1344 !! ** Method : Read the namthd_pnd namelist and check the melt pond 1344 !! ** Method : Read the namthd_pnd namelist and check the melt pond 1345 1345 !! parameter values called at the first timestep (nit000) 1346 1346 !! 1347 !! ** input : Namelist namthd_pnd 1347 !! ** input : Namelist namthd_pnd 1348 1348 !!------------------------------------------------------------------- 1349 1349 INTEGER :: ios, ioptio ! Local integer … … 1389 1389 ! 1390 1390 SELECT CASE( nice_pnd ) 1391 CASE( np_pndNO ) 1391 CASE( np_pndNO ) 1392 1392 IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 1393 1393 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF 1394 CASE( np_pndCST ) 1394 CASE( np_pndCST ) 1395 1395 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF 1396 1396 END SELECT 1397 1397 ! 1398 1398 END SUBROUTINE ice_thd_pnd_init 1399 1399 1400 1400 #else 1401 1401 !!---------------------------------------------------------------------- 1402 1402 !! Default option Empty module NO SI3 sea-ice model 1403 1403 !!---------------------------------------------------------------------- 1404 #endif 1404 #endif 1405 1405 1406 1406 !!====================================================================== 1407 END MODULE icethd_pnd 1407 END MODULE icethd_pnd
Note: See TracChangeset
for help on using the changeset viewer.
