- Timestamp:
- 2020-06-09T18:45:11+02:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE
- Files:
-
- 19 edited
-
ice.F90 (modified) (6 diffs)
-
ice1d.F90 (modified) (4 diffs)
-
icealb.F90 (modified) (1 diff)
-
icedyn.F90 (modified) (1 diff)
-
icedyn_adv.F90 (modified) (1 diff)
-
icedyn_adv_pra.F90 (modified) (13 diffs)
-
icedyn_adv_umx.F90 (modified) (5 diffs)
-
icedyn_rdgrft.F90 (modified) (7 diffs)
-
iceistate.F90 (modified) (8 diffs)
-
iceitd.F90 (modified) (4 diffs)
-
icerst.F90 (modified) (4 diffs)
-
icesbc.F90 (modified) (1 diff)
-
icestp.F90 (modified) (1 diff)
-
icethd.F90 (modified) (2 diffs)
-
icethd_pnd.F90 (modified) (7 diffs)
-
icethd_zdf_bl99.F90 (modified) (4 diffs)
-
iceupdate.F90 (modified) (2 diffs)
-
icevar.F90 (modified) (6 diffs)
-
icewri.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice.F90
r11715 r13080 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 !!-------------|-------------|---------------------------------|-------| … … 195 196 REAL(wp), PUBLIC :: rn_hpnd !: prescribed pond depth (0<rn_hpnd<1) 196 197 LOGICAL , PUBLIC :: ln_pnd_alb !: melt ponds affect albedo 198 REAL(wp), PUBLIC :: rn_pnd_min !: Minimum ice fraction that contributes to melt ponds 199 REAL(wp), PUBLIC :: rn_pnd_max !: Maximum ice fraction that contributes to melt ponds 200 LOGICAL, PUBLIC :: ln_pnd_totfrac !: Use total ice fraction instead of category ice fraction 201 !: when determining ice fraction that contributes to melt ponds 202 LOGICAL, PUBLIC :: ln_pnd_overflow !: Allow ponds to overflow and have vertical flushing 203 LOGICAL, PUBLIC :: ln_pnd_lids !: Melt ponds can have frozen lids 204 LOGICAL, PUBLIC :: ln_use_pndmass !: Use melt pond mass flux diagnostic, passing it to the ocean for emp 197 205 198 206 ! !!** ice-diagnostics namelist (namdia) ** … … 295 303 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_i !: Ice thickness (m) 296 304 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i !: Ice fractional areas (concentration) 305 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time 297 306 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_i !: Ice volume per unit area (m) 298 307 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s !: Snow volume per unit area (m) … … 331 340 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip !: melt pond volume per grid cell area [m] 332 341 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_frac !: melt pond fraction (a_ip/a_i) 342 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_eff !: melt pond effective fraction (not covered up by lid) (a_ip/a_i) 333 343 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond depth [m] 344 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: lh_ip !: melt pond lid thickness [m] 334 345 335 346 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond concentration … … 435 446 436 447 ii = ii + 1 448 ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) ) 449 450 ii = ii + 1 437 451 ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , & 438 452 & vt_i (jpi,jpj) , vt_s (jpi,jpj) , st_i(jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) , & … … 448 462 449 463 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) ) 464 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), & 465 & lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl) , STAT = ierr(ii) ) 451 466 452 467 ii = ii + 1 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice1d.F90
r11715 r13080 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(:) :: a_ip_eff_1d !: 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: lh_ip_1d !: Ice pond lid thickness [m] 130 132 131 133 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 157 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ip_2d 158 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_ip_2d 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: lh_ip_2d 159 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_su_2d 160 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_2d … … 208 211 & dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) , & 209 212 & 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) , 211 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , 213 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , lh_ip_1d(jpij) , & 214 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij) , & 212 215 & sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d (jpij) , STAT=ierr(ii) ) 213 216 ! … … 226 229 ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) , & 227 230 & 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) , 231 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl), & 229 232 & STAT=ierr(ii) ) 230 233 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icealb.F90
r11715 r13080 96 96 LOGICAL , INTENT(in ) :: ld_pnd_alb ! effect of melt ponds on albedo 97 97 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area) 98 ! This is the effective fraction not covered up by a pond lid 98 99 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_pnd ! melt pond depth 99 100 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_cs ! albedo of ice under clear sky -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icedyn_adv.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icedyn_adv_pra.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icedyn_adv_umx.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icedyn_rdgrft.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/iceistate.F90
r11715 r13080 41 41 ! !! ** namelist (namini) ** 42 42 LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not 43 LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file 43 INTEGER, PUBLIC :: nn_iceini_file ! Ice initialization: 44 ! 0 = Initialise sea ice based on SSTs 45 ! 1 = Initialise sea ice from single category netcdf file 46 ! 2 = Initialise sea ice from multi category restart file 44 47 REAL(wp) :: rn_thres_sst 45 48 REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n … … 48 51 REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s 49 52 ! 50 ! ! if ln_iceini_file = T53 ! ! if nn_iceini_file = 1 51 54 INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read 52 55 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) … … 155 158 a_ip (:,:,:) = 0._wp 156 159 v_ip (:,:,:) = 0._wp 160 lh_ip (:,:,:) = 0._wp 157 161 a_ip_frac(:,:,:) = 0._wp 162 a_ip_eff (:,:,:) = 0._wp 158 163 h_ip (:,:,:) = 0._wp 159 164 ! … … 167 172 IF( ln_iceini ) THEN 168 173 ! !---------------! 169 IF( ln_iceini_file)THEN ! Read a file !174 IF( nn_iceini_file == 1 )THEN ! Read a file ! 170 175 ! !---------------! 171 176 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp … … 362 367 a_ip_frac(:,:,:) = 0._wp 363 368 END WHERE 369 a_ip_eff(:,:,:) = a_ip_frac(:,:,:) 364 370 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 365 371 … … 466 472 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 467 473 ! 468 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &474 NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, & 469 475 & rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 470 476 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & … … 493 499 WRITE(numout,*) ' Namelist namini:' 494 500 WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini 495 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file501 WRITE(numout,*) ' ice initialization from a netcdf file nn_iceini_file = ', nn_iceini_file 496 502 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 497 IF( ln_iceini .AND. .NOT.ln_iceini_file) THEN503 IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN 498 504 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 499 505 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s … … 508 514 ENDIF 509 515 ! 510 IF( ln_iceini_file) THEN ! Ice initialization using input file516 IF( nn_iceini_file == 1 ) THEN ! Ice initialization using input file 511 517 ! 512 518 ! set si structure -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceitd.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icerst.F90
r11715 r13080 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 … … 171 172 INTEGER :: jk 172 173 LOGICAL :: llok 173 INTEGER :: id0, id1, id2, id3, id4 ! local integer174 INTEGER :: id0, id1, id2, id3, id4, id5 ! local integer 174 175 CHARACTER(len=25) :: znam 175 176 CHARACTER(len=2) :: zchar, zchar1 … … 250 251 v_ip(:,:,:) = 0._wp 251 252 ENDIF 253 a_ip_eff(:,:,:) = a_ip(:,:,:) 254 ! melt pond lids 255 id5 = iom_varid( numrir, 'lh_ip' , ldstop = .FALSE. ) 256 IF( id5 > 0 ) THEN 257 CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 258 ELSE 259 lh_ip(:,:,:) = 0._wp 260 ENDIF 252 261 ! fields needed for Met Office (Jules) coupling 253 262 IF( ln_cpl ) THEN … … 274 283 CALL ice_istate( nit000 ) 275 284 ! 276 IF( .NOT.ln_iceini .OR. .NOT.ln_iceini_file) &277 & CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and ln_iceini_file=T')285 IF( .NOT.ln_iceini .OR. nn_iceini_file == 0 ) & 286 & CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and nn_iceini_file=0') 278 287 ! 279 288 ENDIF -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icesbc.F90
r11715 r13080 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_eff, h_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_fix_cpl_retain_melt/src/ICE/icestp.F90
r11715 r13080 252 252 ! 253 253 ! ! Initial sea-ice state 254 IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst 255 CALL ice_istate_init 256 CALL ice_istate( nit000 ) 257 ELSE ! start from a restart file 258 CALL ice_rst_read 254 255 CALL ice_istate_init 256 IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 257 CALL ice_rst_read ! start from a restart file 258 ELSE 259 CALL ice_istate( nit000 ) ! start from rest: sea-ice deduced from sst 259 260 ENDIF 260 261 CALL ice_var_glo2eqv -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd.F90
r11715 r13080 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), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 358 CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 357 359 ! 358 360 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 461 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 462 464 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 466 CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 463 467 ! 464 468 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_pnd.F90
r11715 r13080 38 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 39 39 40 REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 41 40 42 !! * Substitutions 41 43 # include "vectopt_loop_substitute.h90" … … 89 91 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 90 92 a_ip_frac_1d(ji) = rn_apnd 93 a_ip_eff_1d(ji) = rn_apnd 91 94 h_ip_1d(ji) = rn_hpnd 92 95 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 93 96 ELSE 94 97 a_ip_frac_1d(ji) = 0._wp 98 a_ip_eff_1d(ji) = 0._wp 95 99 h_ip_1d(ji) = 0._wp 96 100 a_ip_1d(ji) = 0._wp … … 129 133 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 134 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 ! 135 REAL(wp), PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 136 REAL(wp), PARAMETER :: pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 137 ! 138 REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 132 139 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 133 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 140 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) 134 141 REAL(wp) :: z1_Tp ! inverse reference temperature 135 142 REAL(wp) :: z1_rhow ! inverse freshwater density 136 143 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 144 REAL(wp) :: z1_rhoi ! inverse ice density 137 145 REAL(wp) :: zfac, zdum 138 ! 139 INTEGER :: ji ! loop indices 146 REAL(wp) :: t_grad ! Temperature deficit for refreezing 147 REAL(wp) :: omega_dt ! Time independent accumulated variables used for freezing 148 REAL(wp) :: lh_ip_end ! Lid thickness at end of timestep (temporary variable) 149 REAL(wp) :: zdh_frz ! Amount of melt pond that freezes (m) 150 REAL(wp) :: dh_ip_over ! Pond thickness change due to leaking 151 REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean 152 REAL(wp) :: h_gravity_head ! Height above sea level of the top of the melt pond 153 REAL(wp) :: h_percolation ! Distance between the base of the melt pond and the base of the sea ice 154 REAL(wp) :: Sbr ! Brine salinity 155 REAL(wp), DIMENSION(nlay_i) :: phi ! liquid fraction 156 REAL(wp) :: perm ! Permeability of the sea ice 157 REAL(wp) :: za_ip ! Temporary pond fraction 158 REAL(wp) :: max_h_diff_ts ! Maximum meltpond depth change due to leaking or overflow (m per ts) 159 REAL(wp) :: lfrac_pnd ! The fraction of the meltpond exposed (not inder a frozen lid) 160 161 ! 162 INTEGER :: ji, jk ! loop indices 140 163 !!------------------------------------------------------------------- 141 164 z1_rhow = 1._wp / rhow 142 165 z1_zpnd_aspect = 1._wp / zpnd_aspect 143 166 z1_Tp = 1._wp / zTp 167 z1_rhoi = 1._wp / rhoi 168 169 ! Define time-independent field for use in refreezing 170 omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 144 171 145 172 DO ji = 1, npti 146 ! !--------------------------------! 147 IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! 148 ! !--------------------------------! 149 !--- Remove ponds on thin ice 173 174 ! !----------------------------------------------------! 175 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 176 ! !----------------------------------------------------! 177 !--- Remove ponds on thin ice or tiny ice fractions 150 178 a_ip_1d(ji) = 0._wp 151 179 a_ip_frac_1d(ji) = 0._wp 152 180 h_ip_1d(ji) = 0._wp 181 lh_ip_1d(ji) = 0._wp 153 182 ! !--------------------------------! 154 183 ELSE ! Case ice thickness >= rn_himin ! 155 184 ! !--------------------------------! 156 185 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 157 ! 158 ! 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) 160 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 161 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 186 187 ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 188 za_ip = a_ip_1d(ji) 189 IF ( za_ip < epsi06 ) za_ip = epsi06 190 ! 191 ! available meltwater for melt ponding [m, >0] and fraction 192 tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 193 IF ( ln_pnd_totfrac ) THEN 194 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 195 ELSE 196 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji) ! Use category ice fraction 197 ENDIF 198 zdv_mlt = zfr_mlt * tot_mlt 162 199 ! 163 200 !--- 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 ) 201 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 202 ! 203 !--- Lid shrinking. ---! 204 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 166 205 ! 167 206 ! melt pond mass flux (<0) 168 IF( zdv_mlt > 0._wp ) THEN169 zfac = z fr_mlt * zdv_mlt * rhow * r1_rdtice207 IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 208 zfac = zdv_mlt * rhow * r1_rdtice 170 209 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 171 210 ! … … 177 216 ! 178 217 !--- 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 ) 218 IF ( ln_pnd_lids ) THEN 219 220 ! Code to use if using melt pond lids 221 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 222 t_grad = (zTp+rt0) - t_su_1d(ji) 223 224 ! The following equation is a rearranged form of: 225 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 226 ! where: lid_thickness_start = lh_ip_1d(ji) 227 ! lid_thickness_end = lh_ip_end 228 ! omega_dt is a bunch of terms in the equation that do not change 229 ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 230 ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). 231 232 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 233 zdh_frz = lh_ip_end - lh_ip_1d(ji) 234 235 ! Pond shrinking 236 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 237 238 ! Lid growing 239 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 240 ELSE 241 zdh_frz = 0._wp 242 END IF 243 244 ELSE 245 ! Code to use if not using melt pond lids 246 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 247 ENDIF 248 249 ! 250 ! Make sure pond volume or lid thickness has not gone negative 251 IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 252 IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 180 253 ! 181 254 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac … … 184 257 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 185 258 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 259 260 !--- Pond overflow and vertical flushing ---! 261 IF ( ln_pnd_overflow ) THEN 262 263 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 264 IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 265 h_ip_1d(ji) = zpnd_aspect * zfr_mlt 266 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 267 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 268 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 269 ENDIF 270 271 ! If pond depth exceeds half the ice thickness then reduce the pond volume 272 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 273 h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 274 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 275 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 276 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 277 ENDIF 278 279 !-- Vertical flushing of pond water --! 280 ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth. 281 ! This assumes the pond is sitting on top of the ice. 282 h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 283 284 ! The depth through which water percolates is the distance under the melt pond 285 h_percolation = h_i_1d(ji) - h_ip_1d(ji) 286 287 ! Calculate the permeability of the ice (Assur 1958) 288 DO jk = 1, nlay_i 289 Sbr = - 1.2_wp & 290 - 21.8_wp * (t_i_1d(ji,jk)-rt0) & 291 - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 & 292 - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3 293 phi(jk) = sz_i_1d(ji,jk)/Sbr 294 END DO 295 perm = 3.0e-08_wp * (minval(phi))**3 296 297 ! Do the drainage using Darcy's law 298 IF ( perm > 0._wp ) THEN 299 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number 300 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative 301 302 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 303 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 304 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 305 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 306 ENDIF 307 ENDIF 308 309 ! If lid thickness is ten times greater than pond thickness then remove pond 310 IF ( ln_pnd_lids ) THEN 311 IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 312 a_ip_1d(ji) = 0._wp 313 a_ip_frac_1d(ji) = 0._wp 314 h_ip_1d(ji) = 0._wp 315 lh_ip_1d(ji) = 0._wp 316 v_ip_1d(ji) = 0._wp 317 ENDIF 318 ENDIF 319 320 ! If any of the previous changes has removed all the ice thickness then remove ice area. 321 IF ( h_i_1d(ji) == 0._wp ) THEN 322 a_i_1d(ji) = 0._wp 323 h_s_1d(ji) = 0._wp 324 ENDIF 325 186 326 ! 187 327 ENDIF 328 329 ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid) 330 IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN ! Pond lid is very thin. Expose all the pond. 331 lfrac_pnd = 1.0 332 ELSE 333 IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN ! Pond lid is very thick. Cover all the pond up with ice and nosw. 334 lfrac_pnd = 0.0 335 ELSE ! Pond lid is of intermediate thickness. Expose part of the melt pond. 336 lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 337 ENDIF 338 ENDIF 339 340 ! Calculate the melt pond effective area 341 a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd 342 188 343 END DO 189 344 ! … … 205 360 INTEGER :: ios, ioptio ! Local integer 206 361 !! 207 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 362 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb, & 363 rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac, & 364 ln_use_pndmass 208 365 !!------------------------------------------------------------------- 209 366 ! … … 223 380 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 224 381 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 382 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_pnd_min = ', rn_pnd_min 383 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_pnd_max = ', rn_pnd_max 384 WRITE(numout,*) ' Use total ice fraction instead of category ice fraction ln_pnd_totfrac = ',ln_pnd_totfrac 385 WRITE(numout,*) ' Allow ponds to overflow and have vertical flushing ln_pnd_overflow = ',ln_pnd_overflow 386 WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ',ln_pnd_lids 387 WRITE(numout,*) ' Use melt pond mass flux diagnostic, passing to ocean ln_use_pndmass = ',ln_use_pndmass 225 388 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 226 389 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_zdf_bl99.F90
r11715 r13080 31 31 !!---------------------------------------------------------------------- 32 32 !! NEMO/ICE 4.0 , NEMO Consortium (2018) 33 !! $Id $33 !! $Id: icethd_zdf_bl99.F90 10926 2019-05-03 12:32:10Z clem $ 34 34 !! Software governed by the CeCILL license (see ./LICENSE) 35 35 !!---------------------------------------------------------------------- … … 769 769 ! 770 770 ! --- calculate conduction fluxes (positive downward) 771 771 ! bottom ice conduction flux 772 772 DO ji = 1, npti 773 ! ! surface ice conduction flux 774 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 775 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 776 ! ! bottom ice conduction flux 777 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 773 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 778 774 END DO 779 775 ! surface ice conduction flux 776 IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 777 ! 778 DO ji = 1, npti 779 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 780 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 781 END DO 782 ! 783 ELSEIF( k_cnd == np_cnd_ON ) THEN 784 ! 785 DO ji = 1, npti 786 qcn_ice_top_1d(ji) = qcn_ice_1d(ji) 787 END DO 788 ! 789 ENDIF 790 ! surface ice temperature 791 IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN 792 ! 793 DO ji = 1, npti 794 t_su_1d(ji) = ( qcn_ice_top_1d(ji) & ! calculate surface temperature 795 & + isnow(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 796 & + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) & 797 & ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 798 t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp ) ! cap t_su 799 END DO 800 ! 801 ENDIF 780 802 ! 781 803 ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! … … 785 807 DO ji = 1, npti 786 808 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 787 END DO788 !789 ELSEIF( k_cnd == np_cnd_ON ) THEN790 !791 DO ji = 1, npti792 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)793 809 END DO 794 810 ! … … 856 872 t_i_1d (1:npti,:) = ztiold (1:npti,:) 857 873 qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti) 858 859 !!clem860 ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input861 !clem862 874 ENDIF 863 875 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceupdate.F90
r11715 r13080 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_eff, h_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(:,:,:) … … 279 279 IF( iom_use('hfxcndbot' ) ) CALL iom_put( 'hfxcndbot' , SUM( qcn_ice_bot * a_i_b, dim=3 ) ) ! Bottom conduction flux 280 280 IF( iom_use('hfxcndtop' ) ) CALL iom_put( 'hfxcndtop' , SUM( qcn_ice_top * a_i_b, dim=3 ) ) ! Surface conduction flux 281 IF( iom_use('hfxcndcpl' ) ) CALL iom_put( "hfxcndcpl" , SUM( qcn_ice * a_i_b, dim=3 ) ) ! Conduction flux we are giving it 281 282 282 283 ! controls -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icevar.F90
r11715 r13080 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_fix_cpl_retain_melt/src/ICE/icewri.F90
r11715 r13080 149 149 150 150 ! --- category-dependent fields --- ! 151 IF( iom_use('icelhpnd_cat') ) CALL iom_put( 'icelhpnd_cat', lh_ip * zmsk00l ) ! melt pond lid thickness for categories 151 152 IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0% 152 153 IF( iom_use('iceconc_cat' ) ) CALL iom_put( 'iceconc_cat' , a_i * zmsk00l ) ! area for categories … … 165 166 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac for categories 166 167 IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories 168 169 IF( iom_use('icevol_cat' ) ) CALL iom_put( "icevol_cat" , v_i * zmsk00l ) ! volume for categories 170 IF( iom_use('iceaepnd_cat') ) CALL iom_put( 'iceaepnd_cat', a_ip_eff * zmsk00l ) ! melt pond effective frac for categories 167 171 168 172 !------------------
Note: See TracChangeset
for help on using the changeset viewer.
