Changeset 13442 for NEMO/branches
- Timestamp:
- 2020-08-27T17:13:15+02:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice.F90
r13080 r13442 70 70 !! a_ip | - | Ice pond concentration | | 71 71 !! v_ip | - | Ice pond volume per unit area| m | 72 !! s_ip ! s_ip_1d ! Ice pond salinity ! pss ! 72 73 !! lh_ip ! lh_ip_1d ! Ice pond lid thickness ! m ! 73 74 !! | … … 271 272 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_res !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1] 272 273 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_sub !: salt flux due to ice sublimation [pss.kg.m-2.s-1 => g.m-2.s-1] 274 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_pnd !: salt flux due to melt pond-ocean exchange [pss.kg.m-2.s-1 => g.m-2.s-1] 273 275 274 276 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_bog !: total heat flux causing bottom ice growth [W.m-2] … … 343 345 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond depth [m] 344 346 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: lh_ip !: melt pond lid thickness [m] 347 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: s_ip !: melt pond salinity [pss] 345 348 346 349 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond concentration … … 463 466 ii = ii + 1 464 467 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) )468 & lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl), s_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 466 469 467 470 ii = ii + 1 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/ice1d.F90
r13080 r13442 92 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_lam_1d 93 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_dyn_1d 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_pnd_1d 94 95 95 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d … … 115 116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bog !: Ice bottom accretion [m] 116 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_sub !: Ice surface sublimation [m] 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_mlt !: Snow melt [m] 118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_to_pnd !: Ice surface melt to ponds [m] 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_mlt !: Total snow melt [m] 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_to_pnd !: Snow melt going to ponds [m] 118 121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_snowice !: Snow ice formation [m of ice] 119 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_1d !: Ice bulk salinity [ppt] 120 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_pndfrz !: Salinity of new ice after pond freezing 121 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_i_1d !: 122 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_s_1d !: … … 126 130 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_1d !: 127 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_ip_1d !: 132 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_ip_1d !: Ice pond salinity 128 133 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ip_1d !: 129 134 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_frac_1d !: 130 135 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_eff_1d !: 131 136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: lh_ip_1d !: Ice pond lid thickness [m] 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frz_ip !: Heat flux from freezing melt ponds 132 138 133 139 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 159 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ip_2d 160 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_ip_2d 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: s_ip_2d 161 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: lh_ip_2d 162 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_su_2d … … 210 217 & h_i_1d (jpij) , h_ib_1d (jpij) , h_s_1d (jpij) , & 211 218 & dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) , & 212 & dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new (jpij) , &219 & dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new (jpij) , s_i_pndfrz(jpij) ,& 213 220 & 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) ,&221 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij), s_ip_1d(jpij) , frz_ip(jpij) , & 215 222 & sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d (jpij) , STAT=ierr(ii) ) 216 223 ! … … 230 237 & v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) , & 231 238 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl), & 239 & s_ip_2d(jpij,jpl) , & 232 240 & STAT=ierr(ii) ) 233 241 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icectl.F90
r11715 r13442 94 94 ! salt flux 95 95 pdiag_fs = glob_sum( 'icectl', & 96 & ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + &96 & ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_pnd + & 97 97 & sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 98 98 ! heat flux … … 112 112 ! -- salt diag -- ! 113 113 zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice & 114 & + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + 114 & + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_pnd + & 115 115 & sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 116 116 & - pdiag_fs … … 242 242 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 243 243 ! salt flux 244 pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam 244 pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam + sfx_pnd 245 245 ! heat flux 246 246 pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & … … 258 258 ! -- salt diag -- ! 259 259 zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice & 260 & + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &260 & + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam + sfx_pnd ) & 261 261 & - pdiag_fs 262 262 IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel ) ll_stop_s = .TRUE. -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv.F90
r13080 r13442 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, lh_ip, e_s, e_i )86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, sv_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, lh_ip, e_s, e_i )91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, sv_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
r13080 r13442 45 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume 46 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxlhp, sylhp, sxxlhp, syylhp, sxylhp ! melt pond lid thickness 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsalp, sysalp, sxxsalp, syysalp, sxysalp ! melt pond salinity 47 48 48 49 !! * Substitutions … … 56 57 57 58 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, & 58 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, p e_s, pe_i )59 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 59 60 !!---------------------------------------------------------------------- 60 61 !! ** routine ice_dyn_adv_pra ** … … 80 81 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 81 82 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 83 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_ip ! melt pond salt content 82 84 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 83 85 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 132 134 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 133 135 z0lhp(:,:,jl) = plh_ip(:,:,jl) * e1e2t(:,:) ! Melt pond lid thickness 136 z0smp(:,:,jl) = psv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond salt content 134 137 ENDIF 135 138 END DO … … 172 175 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) !--- melt pond lid thickness 173 176 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 177 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) !--- melt pond salinity 178 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) 174 179 ENDIF 175 180 END DO … … 209 214 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) !--- melt pond lid thickness 210 215 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 216 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) !--- melt pond salinity 217 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smp, sxsalp, sxxsalp, sysalp, syysalp, sxysalp) 211 218 ENDIF 212 219 END DO … … 233 240 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 234 241 plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 242 psv_ip(:,:,jl) = z0smp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 235 243 ENDIF 236 244 END DO … … 239 247 ! Remove negative values (conservation is ensured) 240 248 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 241 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, p e_s, pe_i )249 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 242 250 ! 243 251 ! --- Ensure snow load is not too big --- ! … … 662 670 & sxvp(jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 663 671 & sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl), & 672 & sxsalp(jpi,jpj,jpl) , sysalp (jpi,jpj,jpl), sxxsalp (jpi,jpj,jpl), syysalp (jpi,jpj,jpl), sxysalp (jpi,jpj,jpl), & 664 673 ! 665 674 & sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & … … 780 789 CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp ) 781 790 CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp ) 791 ! ! melt pond salinity 792 CALL iom_get( numrir, jpdom_autoglo, 'sxsalp' , sxsalp ) 793 CALL iom_get( numrir, jpdom_autoglo, 'sysalp' , sysalp ) 794 CALL iom_get( numrir, jpdom_autoglo, 'sxxsalp', sxxsalp ) 795 CALL iom_get( numrir, jpdom_autoglo, 'syysalp', syysalp ) 796 CALL iom_get( numrir, jpdom_autoglo, 'sxysalp', sxysalp ) 782 797 ENDIF 783 798 ! … … 798 813 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 799 814 sxlhp = 0._wp ; sylhp = 0._wp ; sxxlhp = 0._wp ; syylhp = 0._wp ; sxylhp = 0._wp ! melt pond lid thickness 815 sxsalp = 0._wp ; sysalp = 0._wp ; sxxsalp = 0._wp ; syysalp = 0._wp ; sxysalp = 0._wp ! melt pond salinity 800 816 ENDIF 801 817 ENDIF … … 884 900 CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp ) 885 901 CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp ) 902 ! ! melt pond salinity 903 CALL iom_rstput( iter, nitrst, numriw, 'sxsalp' , sxsalp ) 904 CALL iom_rstput( iter, nitrst, numriw, 'sysalp' , sysalp ) 905 CALL iom_rstput( iter, nitrst, numriw, 'sxxsalp', sxxsalp ) 906 CALL iom_rstput( iter, nitrst, numriw, 'syysalp', syysalp ) 907 CALL iom_rstput( iter, nitrst, numriw, 'sxysalp', sxysalp ) 886 908 ENDIF 887 909 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_adv_umx.F90
r13080 r13442 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, plh_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, psv_ip, pe_s, pe_i ) 63 63 !!---------------------------------------------------------------------- 64 64 !! *** ROUTINE ice_dyn_adv_umx *** … … 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 87 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 88 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_ip ! melt pond salt content 88 89 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 89 90 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 340 341 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 341 342 & zhvar, plh_ip, zua_ups, zva_ups ) 343 ! salt content 344 zamsk = 0._wp 345 zhvar(:,:,:) = psv_ip(:,:,:) * z1_aip(:,:,:) 346 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 347 & zhvar, psv_ip, zua_ups, zva_ups ) 342 348 343 349 ENDIF … … 357 363 ! Remove negative values (conservation is ensured) 358 364 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 359 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, p e_s, pe_i )365 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 360 366 ! 361 367 ! Make sure ice thickness is not too big -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icedyn_rdgrft.F90
r13080 r13442 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, lhprdg ! area etc of new ridges506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft ! area etc of rafted ice505 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg, svprdg ! area etc of new ridges 506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft, svprft ! area etc of rafted ice 507 507 ! 508 508 REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges … … 579 579 vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 580 580 lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg 581 svprdg(ji) = sv_ip_2d(ji,jl1) * afrdg 581 582 aprft1 = a_ip_2d(ji,jl1) * afrft 582 583 aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 583 584 vprft (ji) = v_ip_2d(ji,jl1) * afrft 584 585 lhprft(ji) = lh_ip_2d(ji,jl1) * afrft 586 svprft(ji) = sv_ip_2d(ji,jl1) * afrft 585 587 ENDIF 586 588 … … 613 615 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 614 616 lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji) 617 sv_ip_2d(ji,jl1) = sv_ip_2d(ji,jl1) - svprdg(ji) - svprft(ji) 615 618 ENDIF 616 619 ENDIF … … 711 714 lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + ( lhprdg (ji) * rn_fpndrdg * fvol (ji) & 712 715 & + lhprft (ji) * rn_fpndrft * zswitch(ji) ) 716 sv_ip_2d (ji,jl2) = sv_ip_2d(ji,jl2) + ( svprdg (ji) * rn_fpndrdg * fvol (ji) & 717 & + svprft (ji) * rn_fpndrft * zswitch(ji) ) 713 718 ENDIF 714 719 … … 741 746 !---------------- 742 747 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 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 )748 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, ,sv_ip_2d, ze_s_2d, ze_i_2d ) 744 749 ! 745 750 END SUBROUTINE rdgrft_shift … … 854 859 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 855 860 CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 861 CALL tab_3d_2d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) ) 856 862 DO jl = 1, jpl 857 863 DO jk = 1, nlay_s … … 868 874 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 869 875 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 876 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd (:,:) ) 870 877 ! 871 878 ! !---------------------! … … 881 888 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 882 889 CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 890 CALL tab_2d_3d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) ) 883 891 DO jl = 1, jpl 884 892 DO jk = 1, nlay_s … … 895 903 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 896 904 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 905 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd (:,:) ) 897 906 ! 898 907 END SELECT -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceistate.F90
r13080 r13442 162 162 a_ip_eff (:,:,:) = 0._wp 163 163 h_ip (:,:,:) = 0._wp 164 s_ip (:,:,:) = 0._wp 164 165 ! 165 166 ! ice velocities -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceitd.F90
r13080 r13442 412 412 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 413 413 CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 414 CALL tab_3d_2d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip ) 414 415 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 415 416 DO jl = 1, jpl … … 488 489 lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans 489 490 lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans 491 ! 492 ztrans = sv_ip_2d(ji,jl1) * zworka(ji) ! Pond salt content 493 sv_ip_2d(ji,jl1) = sv_ip_2d(ji,jl1) - ztrans 494 sv_ip_2d(ji,jl2) = sv_ip_2d(ji,jl2) + ztrans 490 495 ENDIF 491 496 ! … … 532 537 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 533 538 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 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 )539 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, sv_ip_2d, ze_s_2d, ze_i_2d ) 535 540 536 541 ! at_i must be <= rn_amax … … 561 566 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 562 567 CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 568 CALL tab_2d_3d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip ) 563 569 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 564 570 DO jl = 1, jpl -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icerst.F90
r13080 r13442 133 133 CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip ) 134 134 CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip ) 135 CALL iom_rstput( iter, nitrst, numriw, 'sv_ip', sv_ip ) 135 136 ! Snow enthalpy 136 137 DO jk = 1, nlay_s … … 259 260 lh_ip(:,:,:) = 0._wp 260 261 ENDIF 262 ! melt pond salinity 263 id5 = iom_varid( numrir, 'sv_ip' , ldstop = .FALSE. ) 264 IF( id5 > 0 ) THEN 265 CALL iom_get( numrir, jpdom_autoglo, 'sv_ip', sv_ip) 266 ELSE 267 sv_ip(:,:,:) = 0._wp 268 ENDIF 261 269 ! fields needed for Met Office (Jules) coupling 262 270 IF( ln_cpl ) THEN -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icestp.F90
r13080 r13442 404 404 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 405 405 sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp 406 sfx_pnd(:,:) = 0._wp 406 407 ! 407 408 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd.F90
r13080 r13442 217 217 ! 218 218 s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here) 219 s_i_pndfrz(1:npti) = 0._wp 219 220 dh_i_sum (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm (1:npti) = 0._wp 220 221 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 221 222 dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 223 dh_s_to_pnd(1:npti) = 0._wp ; dh_i_to_pnd(1:npti) = 0._wp 222 224 ! 223 225 CALL ice_thd_zdf ! --- Ice-Snow temperature --- ! … … 357 359 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 358 360 CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 361 CALL tab_2d_1d( npti, nptidx(1:npti), s_ip_1d (1:npti), s_ip (:,:,kl) ) 359 362 ! 360 363 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 396 399 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub ) 397 400 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam ) 401 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd ) 398 402 ! 399 403 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd ) … … 465 469 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 466 470 CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d (1:npti), lh_ip (:,:,kl) ) 471 CALL tab_1d_2d( npti, nptidx(1:npti), s_ip_1d (1:npti), s_ip (:,:,kl) ) 467 472 ! 468 473 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 490 495 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub ) 491 496 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam ) 497 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd ) 492 498 ! 493 499 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd ) -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_dh.F90
r11715 r13442 24 24 USE lib_mpp ! MPP library 25 25 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 USE icethd_pnd, ONLY : a_pnd_avail 26 27 27 28 IMPLICIT NONE … … 98 99 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 99 100 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 101 REAL(wp) :: zdeltah_to_ocn 100 102 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 101 103 … … 112 114 CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile 113 115 END SELECT 116 117 ! Initialise fraction of sea ice available for melt ponding 118 DO ji = 1, npti 119 a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji) 120 END DO 114 121 115 122 ! initialize layer thicknesses and enthalpies … … 241 248 zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 242 249 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 250 zdeltah_to_ocn = zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji)) ! Only some of this melt makes it into the ocean. 251 dh_s_to_pnd(ji) = dh_s_to_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji) ! The rest goes onto melt ponds. 252 243 253 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 244 254 245 255 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice ! heat used to melt snow(W.m-2, >0) 246 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah (ji,jk)* r1_rdtice ! snow melting only = water into the ocean (then without snow precip)256 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah_to_ocn * r1_rdtice ! snow melting only = water into the ocean (then without snow precip) 247 257 248 258 ! updates available heat + thickness … … 345 355 346 356 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 357 zdeltah_to_ocn = zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji)) ! Only some of this melt makes it into the ocean. 358 dh_i_to_pnd(ji) = dh_i_to_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji) ! The rest goes onto melt ponds. 347 359 348 360 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat … … 350 362 dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt 351 363 352 zfmdt = - rhoi * zdeltah (ji,jk)! Recompute mass flux [kg/m2, >0]364 zfmdt = - rhoi * zdeltah_to_ocn ! Recompute mass flux [kg/m2, >0] 353 365 354 366 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 355 367 356 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah (ji,jk)* s_i_1d(ji) * r1_rdtice ! Salt flux >0368 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah_to_ocn * s_i_1d(ji) * r1_rdtice ! Salt flux >0 357 369 ! using s_i_1d and not sz_i_1d(jk) is ok) 358 370 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux [W.m-2], < 0 359 371 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], > 0 360 372 ! 361 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah (ji,jk)* r1_rdtice ! Mass flux373 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah_to_ocn * r1_rdtice ! Mass flux 362 374 363 375 END IF … … 403 415 END DO 404 416 417 ! Melt pond freezing 418 !------------------ 419 ! Melt pond freezing amounts already calculated in ice_thd_pnd_frz 420 421 ! update ice enthalpy, thickness and salinity 422 DO ji = 1, npti 423 s_i_pndfrz(ji) = (h_i_1d(ji) * s_i_1d(ji) + dh_i_from_pnd(ji) + s_ip_1d(ji)) / (h_i_1d(ji) + dh_i_from_pnd(ji)) 424 ztmelts = - rTmlt * sz_i_1d(ji,1) ! Freezing point [C] 425 zEi = rcp * ztmelts ! New ice is at freezing point [K] 426 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_i_from_pnd(ji) * (-zEi * rhoi) ! New enthalpy 427 428 429 END DO 430 405 431 406 432 ! Ice Basal growth … … 540 566 ! -------------------------- 541 567 DO ji = 1, npti 542 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )568 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) + dh_i_from_pnd(ji) ) 543 569 END DO 544 570 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_pnd.F90
r13080 r13442 37 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant pond scheme 38 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 39 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 39 40 40 41 REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp … … 132 133 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum - - - - - 133 134 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 134 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature135 135 REAL(wp), PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 136 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 137 ! 138 REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 139 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 138 REAL(wp) :: zfr_mlt ! Total ice and snow surface melt feeding into poinds 140 139 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) 141 140 REAL(wp) :: z1_Tp ! inverse reference temperature … … 180 179 h_ip_1d(ji) = 0._wp 181 180 lh_ip_1d(ji) = 0._wp 181 s_ip_1d(ji) = 0._wp 182 182 ! !--------------------------------! 183 183 ELSE ! Case ice thickness >= rn_himin ! … … 190 190 ! 191 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 192 zdv_mlt = -( dh_i_to_pnd(ji)*rhoi + dh_s_to_pnd(ji)*rhos ) * z1_rhow * a_i_1d(ji) 199 193 ! 200 194 !--- Pond gowth ---! 195 v_ip_old = v_ip_1d(ji) 201 196 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 202 197 ! 203 198 !--- Lid shrinking. ---! 204 199 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 200 ! 201 ! --- Pond salinity changes ---! 202 If ( v_ip_1d(ji) > epsi20 ) THEN 203 s_ip_1d(ji) = ( s_ip_1d(ji) * v_ip_old + s_i_1d(ji) * dh_i_to_pnd(ji)*rhoi* z1_rhow* a_i_1d(ji) ) / v_ip_1d(ji) 204 ELSE 205 s_ip_1d(ji) = 0._wp 206 END IF 205 207 ! 206 208 ! melt pond mass flux (<0) … … 218 220 IF ( ln_pnd_lids ) THEN 219 221 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 222 ! Code to use if using surface heat fluxes to form lids 223 IF ( ln_pnd_hfx ) THEN 224 225 ! Shrink the pond 226 vol_change = dh_i_from_pnd(ji) * a_i_1d(ji) 227 v_ip_1d(ji) = v_ip_1d(ji) - vol_change 228 229 ! Grow the lid 230 lh_ip_1d(ji) = lh_ip_1d(ji) + vol_change / za_ip 231 240 232 ELSE 241 zdh_frz = 0._wp 242 END IF 233 234 ! Code to use if using melt pond lids 235 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 236 t_grad = (zTp+rt0) - t_su_1d(ji) 237 238 ! The following equation is a rearranged form of: 239 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 240 ! where: lid_thickness_start = lh_ip_1d(ji) 241 ! lid_thickness_end = lh_ip_end 242 ! omega_dt is a bunch of terms in the equation that do not change 243 ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 244 ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). 245 246 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 247 zdh_frz = lh_ip_end - lh_ip_1d(ji) 248 249 ! Pond shrinking 250 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 251 252 ! Lid growing 253 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 254 ELSE 255 zdh_frz = 0._wp 256 END IF 257 258 END IF 243 259 244 260 ELSE … … 249 265 ! 250 266 ! 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 267 IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 252 268 IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 253 269 ! … … 263 279 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 264 280 IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 281 v_ip_old = v_ip_1d(ji) 265 282 h_ip_1d(ji) = zpnd_aspect * zfr_mlt 266 283 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 267 284 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 268 285 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 286 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice ! Mass flux to ocean 287 sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 269 288 ENDIF 270 289 271 290 ! If pond depth exceeds half the ice thickness then reduce the pond volume 272 291 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 292 v_ip_old = v_ip_1d(ji) 273 293 h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 274 294 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 275 295 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 276 296 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 297 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice ! Mass flux to ocean 298 sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 277 299 ENDIF 278 300 … … 297 319 ! Do the drainage using Darcy's law 298 320 IF ( perm > 0._wp ) THEN 321 v_ip_old = v_ip_1d(ji) 299 322 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 323 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative … … 304 327 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 305 328 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 329 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * r1_rdtice ! Mass flux to ocean 330 sfx_pnd_1d(ji) = sfx_pnd_1d(ji) + rhow * ( v_ip_old - v_ip_1d(ji)) * s_ip_1d(ji) * r1_rdtice ! Salinity flux to ocean 317 331 ENDIF 318 332 ENDIF … … 322 336 a_i_1d(ji) = 0._wp 323 337 h_s_1d(ji) = 0._wp 338 ENDIF 339 340 ! If there is no melt pond left then set all pond stats to zero 341 IF ( v_ip_1d(ji) <= 0._wp ) THEN 342 v_ip_1d(ji) = 0._wp 343 a_ip_1d(ji) = 0._wp 344 a_ip_frac_1d(ji) = 0._wp 345 h_ip_1d(ji) = 0._wp 346 lh_ip_1d(ji) = 0._wp 347 s_ip_1d(ji) = 0._wp 324 348 ENDIF 325 349 … … 344 368 ! 345 369 END SUBROUTINE pnd_H12 370 371 SUBROUTINE ice_thd_pnd_frz 372 !!-------------------------------------------------------------------------- 373 !! *** ROUTINE ice_thd_pnd_frz *** 374 !! 375 !! ** Purpose : calculate the amount of energy released from pond freezing, 376 !! alter the surface heat flux and generate new ice. 377 !! 378 !!--------------------------------------------------------------------------- 379 380 REAL(wp), PARAMETER :: t_max_frz = -5.0_wp ! Ice surface temperature below which all surface flux goes into freezing 381 382 INTEGER :: ji ! loop indices 383 REAL(wp) :: vol_pnd ! Volume of melt pond per m2 of ice 384 REAL(wp) :: t_based_frac ! Fraction contribution to freezing based on temperature 385 REAL(wp) :: p_based_frac ! Fraction contribution to freezing based on pond area 386 REAL(wp) :: total_frac ! Total fraction of surface heat used in pond freezing 387 REAL(wp) :: frz_whole_pnd ! How much ebergy is released by completely freezing the melt pong 388 389 DO ji = 1, npti 390 391 top_frz_1d(ji) = 0._wp 392 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! Get the pond volume per m2 of grid box 393 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) ! Get the pond fraction 394 vol_pnd = h_ip_1d(ji) * a_ip_frac_1d(ji) ! Get the pond volume per m2 of ice 395 396 ! Only start freezing ponds when the following criteria have been met: 397 IF ( qml_ice_1d(ji) < 0.01_wp .AND. ! There is virtually no top melt 398 qcn_ice_1d(ji) < 0._wp .AND. ! There is a negative surface heat flux (the surface is cooling down) 399 t_su_1d(ji) < (zTp+rt0) .AND. ! The temperature of the surrounding ice is less then the reference temperature 400 vol_pnd > 0._wp ) THEN ! There is a pond to freeze 401 402 ! Calculate the fraction of heat to use based on how cold the surface temperature is 403 t_based_frac = (t_su_1d(ji) - (zTp+rt0))/t_max_frz 404 IF ( t_based_frac > 1._wp ) t_based_frac = 1._wp 405 406 ! Calculate the fraction of heat to use based on how much pond there is. 407 ! Use a fractional area 0.1 larger as surrounding ice can contribute to pond freezing 408 p_based_frac = a_ip_frac_1d(ji) + 0.1_wp 409 IF ( p_based_frac > 1._wp ) p_based_frac = 1._wp 410 411 ! Combine both fractional contributions 412 total_frac = p_based_frac * t_based_frac 413 414 ! What is the total amount of energy released by completely freezing the melt pond 415 frz_whole_pnd = vol_pnd * rhow * rLfus 416 417 ! How much energy is released from freezing 418 frz_ip(ji) = -1.0 * qcn_ice_1d(ji) * total_frac 419 420 ! Make sure we are not freezng more melt pond than there is 421 IF ( frz_ip(i) * rdt_ice > frz_whole_pnd ) THEN 422 frz_ip(ji) = frz_whole_pnd * r1_rdtice 423 END IF 424 425 ! Alter the surface energy flux to account for some energy from freezing 426 qcn_ice_1d(ji) = qcn_ice_1d(ji) + frz_ip(ji) 427 428 ! How much sea ice growth is caused by freezing melt ponds 429 dh_i_from_pnd(ji) = frz_ip(ji) * rdtice * r1_rLfus * r1_rhoi 430 431 END IF 432 END DO 433 434 END SUBROUTINE ice_thd_pnd_frz 346 435 347 436 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icethd_sal.F90
r11715 r13442 80 80 zs_i_si = ( zs_sni - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice 81 81 zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog (ji) / h_i_1d(ji) ! bottom growth 82 zs_i_pnd = (s_i_pndfrz(ji) - s_i_1d(ji) ) * dh_i_from_pnd (ji) / h_i_1d(ji) ! bottom growth 82 83 ! Update salinity (nb: salt flux already included in icethd_dh) 83 s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si 84 s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si + zs_i_pnd 84 85 ENDIF 85 86 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceupdate.F90
r13080 r13442 165 165 !------------------------------------------ 166 166 sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj) & 167 & + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) 167 & + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) & 168 & + sfx_pnd(ji,jj) 168 169 169 170 ! Mass of snow and ice per unit area … … 209 210 IF( iom_use('sfxres' ) ) CALL iom_put( 'sfxres', sfx_res * 1.e-03 ) ! salt flux from undiagnosed processes 210 211 IF( iom_use('sfxsub' ) ) CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 ) ! salt flux from sublimation 212 IF( iom_use('sfxsub' ) ) CALL iom_put( 'sfxpnd', sfx_pnd * 1.e-03 ) ! salt flux from melt ponds 211 213 212 214 ! --- mass fluxes [kg/m2/s] --- ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/icevar.F90
r13080 r13442 534 534 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 535 535 lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj) 536 sv_ip (ji,jj,jl) = sv_ip (ji,jj,jl) * zswitch(ji,jj) 536 537 ! 537 538 END DO … … 556 557 557 558 558 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, p e_s, pe_i )559 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 559 560 !!------------------------------------------------------------------- 560 561 !! *** ROUTINE ice_var_zapneg *** … … 572 573 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 573 574 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 575 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_ip ! melt pond salt content 574 576 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 575 577 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 639 641 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 640 642 ! but it does not change conservation, so keep it this way is ok 643 WHERE( psv_ip (:,:,:) < 0._wp ) psv_ip (:,:,:) = 0._wp 641 644 WHERE( plh_ip (:,:,:) < 0._wp ) plh_ip (:,:,:) = 0._wp 642 645 ! … … 644 647 645 648 646 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, p e_s, pe_i )649 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, psv_ip, pe_s, pe_i ) 647 650 !!------------------------------------------------------------------- 648 651 !! *** ROUTINE ice_var_roundoff *** … … 658 661 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 659 662 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: plh_ip ! melt pond lid thickness 663 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_ip ! melt pond salt content 660 664 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 661 665 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 673 677 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 678 WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 ) plh_ip(1:npti,:) = 0._wp ! lh_ip must be >= 0 679 WHERE( psv_ip(1:npti,:) < 0._wp .AND. psv_ip(1:npti,:) > -epsi10 ) psv_ip(1:npti,:) = 0._wp ! sv_ip must be >= 0 675 680 ENDIF 676 681 ! … … 792 797 !!------------------------------------------------------------------- 793 798 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 794 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )799 & ptmi, ptms, ptmsu, psmi, patip, phtip, psmip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 795 800 !!------------------------------------------------------------------- 796 801 !! ** Purpose : converting 1-cat ice to 1 ice category … … 798 803 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 799 804 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 800 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds801 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds805 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, psmip ! input ice/snow temp & sal & ponds 806 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ! output ice/snow temp & sal & ponds 802 807 !!------------------------------------------------------------------- 803 808 ! == thickness and concentration == ! … … 813 818 pa_ip(:) = patip(:) 814 819 ph_ip(:) = phtip(:) 820 ps_ip (:) = psmip (:) 815 821 816 822 END SUBROUTINE ice_var_itd_1c1c 817 823 818 824 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 819 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )825 & ptmi, ptms, ptmsu, psmi, patip, phtip, psmip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 820 826 !!------------------------------------------------------------------- 821 827 !! ** Purpose : converting N-cat ice to 1 ice category … … 823 829 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 824 830 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 825 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds826 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds831 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, psmip ! input ice/snow temp & sal & ponds 832 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ! output ice/snow temp & sal & ponds 827 833 ! 828 834 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs … … 862 868 ELSEWHERE ; ph_ip(:) = 0._wp 863 869 END WHERE 870 WHERE( pa_ip(:) /= 0._wp ) ; ps_ip(:) = SUM( psmip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 871 ELSEWHERE ; ps_ip(:) = 0._wp 872 END WHERE 864 873 ! 865 874 DEALLOCATE( z1_ai, z1_vi, z1_vs ) … … 868 877 869 878 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 870 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )879 & ptmi, ptms, ptmsu, psmi, patip, phtip, psmip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 871 880 !!------------------------------------------------------------------- 872 881 !! … … 890 899 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 891 900 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 892 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds893 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds901 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, psmip ! input ice/snow temp & sal & ponds 902 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ! output ice/snow temp & sal & ponds 894 903 ! 895 904 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti … … 1001 1010 ELSEWHERE ; ph_ip(:,jl) = 0._wp 1002 1011 END WHERE 1012 ps_ip (:,jl) = psmip (:) 1003 1013 END DO 1004 1014 DEALLOCATE( zfra ) … … 1007 1017 1008 1018 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1009 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )1019 & ptmi, ptms, ptmsu, psmi, patip, phtip, psmip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ) 1010 1020 !!------------------------------------------------------------------- 1011 1021 !! … … 1038 1048 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1039 1049 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1040 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds1041 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds1050 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, psmip ! input ice/snow temp & sal & ponds 1051 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ps_ip ! output ice/snow temp & sal & ponds 1042 1052 ! 1043 1053 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 … … 1068 1078 pa_ip(:,:) = patip(:,:) 1069 1079 ph_ip(:,:) = phtip(:,:) 1080 ps_ip(:,:) = psmip (:,:) 1070 1081 ! ! ---------------------- ! 1071 1082 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! … … 1073 1084 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 1074 1085 & ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1075 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &1076 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) )1086 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), psmip(:,1), & 1087 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ps_ip(:,:) ) 1077 1088 ! ! ---------------------- ! 1078 1089 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! … … 1080 1091 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 1081 1092 & ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1082 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &1083 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) )1093 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), psmip(:,:), & 1094 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ps_ip(:,1) ) 1084 1095 ! ! ----------------------- ! 1085 1096 ELSE ! input cat /= output cat ! … … 1201 1212 END DO 1202 1213 ! 1203 DEALLOCATE( z1_ai, z1_vi, z1_vs , ztmp)1214 DEALLOCATE( z1_ai, z1_vi, z1_vs ) 1204 1215 ! 1205 1216 ! == ponds == ! … … 1223 1234 END WHERE 1224 1235 END DO 1225 DEALLOCATE( zfra ) 1236 ztmp(:) = SUM( psmip (:,:) * patip(:,:) * phtip(:,:), dim=2 ) / SUM( phtip(:,:) * patip(:,:), dim=2 ) 1237 DO jl = 1, jpl 1238 ps_i (:,jl) = ztmp(:) 1239 END DO 1240 DEALLOCATE( zfra, ztmp ) 1226 1241 ! 1227 1242 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.