Changeset 12401
- Timestamp:
- 2020-02-18T16:21:31+01:00 (5 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src
- Files:
-
- 18 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice.F90
r11715 r12401 70 70 !! a_ip | - | Ice pond concentration | | 71 71 !! v_ip | - | Ice pond volume per unit area| m | 72 !! lh_ip ! lh_ip_1d ! Ice pond lid thickness ! m ! 72 73 !! | 73 74 !!-------------|-------------|---------------------------------|-------| … … 332 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_frac !: melt pond fraction (a_ip/a_i) 333 334 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond depth [m] 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: lh_ip !: melt pond lid thickness [m] 334 336 335 337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond concentration … … 448 450 449 451 ii = ii + 1 450 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )452 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), lh_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 451 453 452 454 ii = ii + 1 -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice1d.F90
r11715 r12401 128 128 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ip_1d !: 129 129 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_frac_1d !: 130 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: lh_ip_1d !: Ice pond lid thickness [m] 130 131 131 132 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 157 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ip_2d 158 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_ip_2d 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: lh_ip_2d 159 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_su_2d 160 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_2d … … 208 210 & dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) , & 209 211 & dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new (jpij) , & 210 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , 212 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , lh_ip_1d(jpij) , & 211 213 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , & 212 214 & sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d (jpij) , STAT=ierr(ii) ) … … 226 228 ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) , & 227 229 & v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) , & 228 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , 230 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl), & 229 231 & STAT=ierr(ii) ) 230 232 -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icealb.F90
r11715 r12401 45 45 CONTAINS 46 46 47 SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, p alb_cs, palb_os )47 SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, plh_pnd, palb_cs, palb_os ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE ice_alb *** … … 97 97 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area) 98 98 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_pnd ! melt pond depth 99 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: plh_pnd ! melt pond lid thickness 99 100 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_cs ! albedo of ice under clear sky 100 101 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_os ! albedo of ice under overcast sky 102 ! 103 REAL(wp), PARAMETER :: pnd_lid_max = 0.015 ! pond lid thickness above which the ponds disappear from the albedo calculation 104 REAL(wp), PARAMETER :: pnd_lid_min = 0.005 ! pond lid thickness below which the full pond area is used in the albedo calculation 105 ! Note: these two variables are mirrored in sbccpl.F90 (maybe put them in one place...) 101 106 ! 102 107 INTEGER :: ji, jj, jl ! dummy loop indices … … 106 111 REAL(wp) :: zalb_ice, zafrac_ice ! bare sea ice albedo & relative ice fraction 107 112 REAL(wp) :: zalb_snw, zafrac_snw ! snow-covered sea ice albedo & relative snow fraction 113 REAL(wp) :: lfrac_pnd ! The fraction of the meltpond exposed (not inder a frozen lid) 108 114 !!--------------------------------------------------------------------- 109 115 ! … … 123 129 zafrac_snw = 0._wp 124 130 IF( ld_pnd_alb ) THEN 125 zafrac_pnd = pafrac_pnd(ji,jj,jl) 131 IF ( plh_pnd(ji,jj,jl) > pnd_lid_max ) THEN 132 lfrac_pnd = 0._wp 133 ELSE 134 IF ( plh_pnd(ji,jj,jl) < pnd_lid_min ) THEN 135 lfrac_pnd = 1._wp 136 ELSE 137 lfrac_pnd = ( plh_pnd(ji,jj,jl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 138 END IF 139 END IF 140 zafrac_pnd = pafrac_pnd(ji,jj,jl) * lfrac_pnd 126 141 ELSE 127 142 zafrac_pnd = 0._wp -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn.F90
r11715 r12401 101 101 ELSEWHERE 102 102 h_ip(:,:,:) = 0._wp 103 lh_ip(:,:,:) = 0._wp 103 104 END WHERE 104 105 ! -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv.F90
r11715 r12401 84 84 ! !-----------------------! 85 85 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 87 87 ! !-----------------------! 88 88 CASE( np_advPRA ) ! PRATHER scheme ! 89 89 ! !-----------------------! 90 90 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, & 91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 92 92 END SELECT 93 93 -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_pra.F90
r11715 r12401 44 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction 45 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxlhp, sylhp, sxxlhp, syylhp, sxylhp ! melt pond lid thickness 46 47 47 48 !! * Substitutions … … 55 56 56 57 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, & 57 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )58 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 58 59 !!---------------------------------------------------------------------- 59 60 !! ** routine ice_dyn_adv_pra ** … … 78 79 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 79 80 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 81 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 80 82 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 81 83 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 88 90 REAL(wp), DIMENSION(jpi,jpj,1) :: z0opw 89 91 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 90 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0lhp 91 93 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 92 94 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei … … 129 131 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 130 132 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 133 z0lhp(:,:,jl) = plh_ip(:,:,jl) * e1e2t(:,:) ! Melt pond lid thickness 131 134 ENDIF 132 135 END DO … … 167 170 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 168 171 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 172 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) !--- melt pond lid thickness 173 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 169 174 ENDIF 170 175 END DO … … 202 207 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 203 208 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 209 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) !--- melt pond lid thickness 210 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 204 211 ENDIF 205 212 END DO … … 225 232 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 226 233 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 234 plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 227 235 ENDIF 228 236 END DO … … 231 239 ! Remove negative values (conservation is ensured) 232 240 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 233 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )241 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 234 242 ! 235 243 ! --- Ensure snow load is not too big --- ! … … 653 661 & sxap(jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , & 654 662 & sxvp(jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 663 & sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl), & 655 664 ! 656 665 & sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & … … 765 774 CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 766 775 CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 776 ! ! melt pond lid thickness 777 CALL iom_get( numrir, jpdom_autoglo, 'sxlhp' , sxlhp ) 778 CALL iom_get( numrir, jpdom_autoglo, 'sylhp' , sylhp ) 779 CALL iom_get( numrir, jpdom_autoglo, 'sxxlhp', sxxlhp ) 780 CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp ) 781 CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp ) 767 782 ENDIF 768 783 ! … … 782 797 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 783 798 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 799 sxlhp = 0._wp ; sylhp = 0._wp ; sxxlhp = 0._wp ; syylhp = 0._wp ; sxylhp = 0._wp ! melt pond lid thickness 784 800 ENDIF 785 801 ENDIF … … 862 878 CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 863 879 CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 880 ! ! melt pond lid thickness 881 CALL iom_rstput( iter, nitrst, numriw, 'sxlhp' , sxlhp ) 882 CALL iom_rstput( iter, nitrst, numriw, 'sylhp' , sylhp ) 883 CALL iom_rstput( iter, nitrst, numriw, 'sxxlhp', sxxlhp ) 884 CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp ) 885 CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp ) 864 886 ENDIF 865 887 ! -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_umx.F90
r11715 r12401 60 60 61 61 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 63 63 !!---------------------------------------------------------------------- 64 64 !! *** ROUTINE ice_dyn_adv_umx *** … … 85 85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 87 88 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 88 89 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 334 335 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 335 336 & zhvar, pv_ip, zua_ups, zva_ups ) 337 ! lid thickness 338 zamsk = 0._wp 339 zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:) 340 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 341 & zhvar, plh_ip, zua_ups, zva_ups ) 342 336 343 ENDIF 337 344 ! … … 350 357 ! Remove negative values (conservation is ensured) 351 358 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 352 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )359 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 353 360 ! 354 361 ! Make sure ice thickness is not too big … … 1519 1526 !! 1520 1527 !! ** Purpose : Thickness correction in case advection scheme creates 1521 !! abnormally t ick ice or snow1528 !! abnormally thick ice or snow 1522 1529 !! 1523 1530 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_rdgrft.F90
r11715 r12401 503 503 REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1 504 504 REAL(wp) :: airft1, oirft1, aprft1 505 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg ! area etc of new ridges506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft ! area etc of rafted ice505 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg ! area etc of new ridges 506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft ! area etc of rafted ice 507 507 ! 508 508 REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges … … 578 578 aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 579 579 vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 580 lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg 580 581 aprft1 = a_ip_2d(ji,jl1) * afrft 581 582 aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 582 583 vprft (ji) = v_ip_2d(ji,jl1) * afrft 584 lhprft(ji) = lh_ip_2d(ji,jl1) * afrft 583 585 ENDIF 584 586 … … 610 612 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 611 613 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 614 lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji) 612 615 ENDIF 613 616 ENDIF … … 706 709 a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & 707 710 & + aprft2(ji) * rn_fpndrft * zswitch(ji) ) 711 lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + ( lhprdg (ji) * rn_fpndrdg * fvol (ji) & 712 & + lhprft (ji) * rn_fpndrft * zswitch(ji) ) 708 713 ENDIF 709 714 … … 736 741 !---------------- 737 742 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 738 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )743 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 739 744 ! 740 745 END SUBROUTINE rdgrft_shift … … 848 853 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 849 854 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 855 CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 850 856 DO jl = 1, jpl 851 857 DO jk = 1, nlay_s … … 874 880 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 875 881 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 882 CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 876 883 DO jl = 1, jpl 877 884 DO jk = 1, nlay_s -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceistate.F90
r11715 r12401 155 155 a_ip (:,:,:) = 0._wp 156 156 v_ip (:,:,:) = 0._wp 157 lh_ip (:,:,:) = 0._wp 157 158 a_ip_frac(:,:,:) = 0._wp 158 159 h_ip (:,:,:) = 0._wp -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceitd.F90
r11715 r12401 411 411 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 412 412 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 413 CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 413 414 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 414 415 DO jl = 1, jpl … … 483 484 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 484 485 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 486 ! 487 ztrans = lh_ip_2d(ji,jl1) * zworka(ji) ! Pond lid thickness 488 lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans 489 lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans 485 490 ENDIF 486 491 ! … … 527 532 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 528 533 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 529 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )534 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 530 535 531 536 ! at_i must be <= rn_amax … … 555 560 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 556 561 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 562 CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 557 563 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 558 564 DO jl = 1, jpl -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icerst.F90
r11715 r12401 132 132 CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip ) 133 133 CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip ) 134 CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip ) 134 135 ! Snow enthalpy 135 136 DO jk = 1, nlay_s … … 245 246 CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip ) 246 247 CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip ) 248 CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 247 249 ELSE ! start from rest 248 250 IF(lwp) WRITE(numout,*) ' ==>> previous run without melt ponds output then set it to zero' 249 251 a_ip(:,:,:) = 0._wp 250 252 v_ip(:,:,:) = 0._wp 253 lh_ip(:,:,:) = 0._wp 251 254 ENDIF 252 255 ! fields needed for Met Office (Jules) coupling -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icesbc.F90
r11715 r12401 132 132 133 133 ! --- cloud-sky and overcast-sky ice albedos --- ! 134 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )134 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, lh_ip, zalb_cs, zalb_os ) 135 135 136 136 ! albedo depends on cloud fraction because of non-linear spectral effects -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd.F90
r11715 r12401 355 355 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 356 356 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 357 CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 357 358 ! 358 359 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 461 462 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 462 463 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 463 465 ! 464 466 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90
r11715 r12401 129 129 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 130 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 REAL(wp), PARAMETER :: freezing_t = 273.0_wp ! temperature below which refreezing occurs 131 132 ! 132 133 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 133 134 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 135 REAL(wp) :: actual_mlt ! actual melt water used to make melt ponds (per m2). 136 ! Need to multiply this by ice area to work on volumes. 134 137 REAL(wp) :: z1_Tp ! inverse reference temperature 135 138 REAL(wp) :: z1_rhow ! inverse freshwater density 136 139 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 137 140 REAL(wp) :: zfac, zdum 141 REAL(wp) :: t_grad ! Temperature deficit for refreezing 142 REAL(wp) :: omega_dt ! Time independent accumulated variables used for freezing 143 REAL(wp) :: lh_ip_end ! Lid thickness at end of timestep (temporary variable) 144 REAL(wp) :: actual_frz ! Amount of melt pond that freezes 138 145 ! 139 146 INTEGER :: ji ! loop indices … … 142 149 z1_zpnd_aspect = 1._wp / zpnd_aspect 143 150 z1_Tp = 1._wp / zTp 151 152 ! Define time-independent field for use in refreezing 153 omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 144 154 145 155 DO ji = 1, npti … … 151 161 a_ip_frac_1d(ji) = 0._wp 152 162 h_ip_1d(ji) = 0._wp 163 lh_ip_1d(ji) = 0._wp 164 165 actual_mlt = 0._wp 166 actual_frz = 0._wp 153 167 ! !--------------------------------! 154 168 ELSE ! Case ice thickness >= rn_himin ! … … 157 171 ! 158 172 ! available meltwater for melt ponding [m, >0] and fraction 159 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)173 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 160 174 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 161 175 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 176 actual_mlt = zfr_mlt * zdv_mlt 162 177 ! 163 178 !--- Pond gowth ---! 164 ! v_ip should never be negative, otherwise code crashes 165 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 166 ! 179 v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji) 180 ! 181 !--- Lid shrinking ---! 182 lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt 183 ! 184 ! 185 !--- Pond contraction (due to refreezing) ---! 186 IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN 187 t_grad = freezing_t - t_su_1d(ji) 188 189 ! The following equation is a rearranged form of: 190 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 191 ! where: lid_thickness_start = lh_ip_1d(ji) 192 ! lid_thickness_end = lh_ip_end 193 ! omega_dt is a bunch of terms in the equation that do not change 194 ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same 195 ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification. 196 197 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 198 actual_frz = lh_ip_end - lh_ip_1d(ji) 199 200 ! Pond shrinking 201 v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji) 202 203 ! Lid growing 204 lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz 205 ELSE 206 actual_frz = 0._wp 207 END IF 208 167 209 ! melt pond mass flux (<0) 168 210 IF( zdv_mlt > 0._wp ) THEN 169 zfac = zfr_mlt * zdv_mlt* rhow * r1_rdtice211 zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice 170 212 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 171 213 ! … … 176 218 ENDIF 177 219 ! 178 !--- Pond contraction (due to refreezing) ---! 179 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 220 ! Make sure pond volume or lid thickness has not gone negative 221 IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 222 IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 180 223 ! 181 224 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac … … 184 227 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 185 228 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 229 230 ! If lid thickness is ten times greater than pond thickness then remove pond 231 IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 232 a_ip_1d(ji) = 0._wp 233 a_ip_frac_1d(ji) = 0._wp 234 h_ip_1d(ji) = 0._wp 235 lh_ip_1d(ji) = 0._wp 236 v_ip_1d(ji) = 0._wp 237 END IF 186 238 ! 187 239 ENDIF -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceupdate.F90
r11715 r12401 185 185 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 186 186 !------------------------------------------------------------------ 187 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos187 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, lh_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 188 188 ! 189 189 alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icevar.F90
r11715 r12401 533 533 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 534 534 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 535 lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj) 535 536 ! 536 537 END DO … … 555 556 556 557 557 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )558 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 558 559 !!------------------------------------------------------------------- 559 560 !! *** ROUTINE ice_var_zapneg *** … … 570 571 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 571 572 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 573 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 572 574 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 573 575 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 637 639 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 638 640 ! but it does not change conservation, so keep it this way is ok 641 WHERE( plh_ip (:,:,:) < 0._wp ) plh_ip (:,:,:) = 0._wp 639 642 ! 640 643 END SUBROUTINE ice_var_zapneg 641 644 642 645 643 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, p e_s, pe_i )646 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 644 647 !!------------------------------------------------------------------- 645 648 !! *** ROUTINE ice_var_roundoff *** … … 654 657 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 655 658 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 659 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 656 660 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 657 661 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 668 672 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 669 673 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 674 WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 ) plh_ip(1:npti,:) = 0._wp ! lh_ip must be >= 0 670 675 ENDIF 671 676 ! -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/BDY/bdyice.F90
r11715 r12401 95 95 ! exchange 3d arrays 96 96 CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. & 97 & , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &97 & , a_ip, 'T', 1., v_ip, 'T', 1., lh_ip,'T', 1., s_i , 'T', 1., t_su, 'T', 1. & 98 98 & , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1. & 99 99 & , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) … … 163 163 a_ip(ji,jj, jl) = ( a_ip(ji,jj, jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond concentration 164 164 h_ip(ji,jj, jl) = ( h_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond depth 165 lh_ip(ji,jj, jl) = ( lh_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice pond lid thickness 165 166 ! 166 167 sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) … … 170 171 a_ip(ji,jj,jl) = 0._wp 171 172 h_ip(ji,jj,jl) = 0._wp 173 lh_ip(ji,jj,jl) = 0._wp 172 174 ENDIF 173 175 ! … … 231 233 a_ip(ji,jj, jl) = a_ip(ib,jb, jl) 232 234 h_ip(ji,jj, jl) = h_ip(ib,jb, jl) 235 lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl) 233 236 ! 234 237 sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) … … 280 283 a_ip(ji,jj, jl) = 0._wp 281 284 v_ip(ji,jj, jl) = 0._wp 285 lh_ip(ji,jj, jl) = 0._wp 282 286 t_su(ji,jj, jl) = rt0 283 287 t_s (ji,jj,:,jl) = rt0 -
NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/SBC/sbccpl.F90
r11715 r12401 2047 2047 INTEGER, INTENT(in) :: kt 2048 2048 ! 2049 REAL(wp), PARAMETER :: pnd_lid_max = 0.015 ! pond lid thickness above which the ponds disappear from the albedo calculation 2050 REAL(wp), PARAMETER :: pnd_lid_min = 0.005 ! pond lid thickness below which the full pond area is used in the albedo calculation 2051 ! Note: these two variables are mirrored in icealb.F90 (maybe put them in one place...) 2052 ! 2049 2053 INTEGER :: ji, jj, jl ! dummy loop indices 2050 2054 INTEGER :: isec, info ! local integer … … 2052 2056 REAL(wp), DIMENSION(jpi,jpj) :: zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 2053 2057 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztmp3, ztmp4 2058 REAL(wp), DIMENSION(jpi,jpj,jpl) :: lfrac_pnd ! The fraction of the meltpond exposed (not inder a frozen lid) 2059 2054 2060 !!---------------------------------------------------------------------- 2055 2061 ! … … 2256 2262 SELECT CASE( sn_snd_mpnd%clcat ) 2257 2263 CASE( 'yes' ) 2258 ztmp3(:,:,1:jpl) = a_ip(:,:,1:jpl) 2259 ztmp4(:,:,1:jpl) = v_ip(:,:,1:jpl) 2264 2265 ! Calculate how much meltpond is exposed (not under a frozen lid) 2266 lfrac_pnd(:,:,1:jpl) = 1.0 2267 WHERE( lh_ip(:,:,1:jpl) > pnd_lid_max ) 2268 lfrac_pnd(:,:,1:jpl) = 0.0 2269 END WHERE 2270 WHERE( lh_ip(:,:,1:jpl) > pnd_lid_min .AND. lh_ip(:,:,1:jpl) <= pnd_lid_max ) 2271 lfrac_pnd(:,:,1:jpl) = ( lh_ip(:,:,1:jpl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 2272 END WHERE 2273 2274 ztmp3(:,:,1:jpl) = a_ip_frac(:,:,1:jpl) * lfrac_pnd(:,:,1:jpl) 2275 ztmp4(:,:,1:jpl) = h_ip(:,:,1:jpl) 2276 2260 2277 CASE( 'no' ) 2261 2278 ztmp3(:,:,:) = 0.0
Note: See TracChangeset
for help on using the changeset viewer.