Changeset 14026
- Timestamp:
- 2020-12-03T09:48:10+01:00 (4 years ago)
- Location:
- NEMO/releases/r4.0/r4.0-HEAD/src/ICE
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90
r13640 r14026 366 366 !!---------------------------------------------------------------------- 367 367 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 368 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip_b, v_il_b !: ponds and lids volumes 368 369 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b !: 369 370 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content … … 392 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vpnd !: pond volume variation [m/s] 394 396 ! 395 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) … … 468 470 & et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) , & 469 471 & sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) , & 470 & om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj) 472 & om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), STAT=ierr(ii) ) 471 473 472 474 ii = ii + 1 … … 486 488 ii = ii + 1 487 489 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), & 490 & v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , & 488 491 & a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 489 492 & STAT=ierr(ii) ) … … 500 503 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 501 504 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 502 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), &505 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 503 506 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 504 507 -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90
r13589 r14026 85 85 !! 86 86 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, & 87 & zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 87 & zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 88 & zdiag_eimin, zdiag_esmin, zdiag_simin 88 89 REAL(wp) :: zvtrp, zetrp 89 90 REAL(wp) :: zarea … … 92 93 IF( icount == 0 ) THEN 93 94 94 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos , dim=3 ) * e1e2t )95 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi 95 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 96 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 96 97 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 97 98 … … 112 113 113 114 ! -- mass diag -- ! 114 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice & 115 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) & 116 & - pdiag_v ) * r1_rdtice & 115 117 & + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + & 116 118 & wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & … … 132 134 133 135 ! -- min/max diag -- ! 134 zdiag_amax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 135 zdiag_vmin = glob_min( 'icectl', v_i ) 136 zdiag_amin = glob_min( 'icectl', a_i ) 137 zdiag_smin = glob_min( 'icectl', sv_i ) 136 zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 137 zdiag_vimin = glob_min( 'icectl', v_i ) 138 zdiag_vsmin = glob_min( 'icectl', v_s ) 139 zdiag_vpmin = glob_min( 'icectl', v_ip ) 140 zdiag_vlmin = glob_min( 'icectl', v_il ) 141 zdiag_aimin = glob_min( 'icectl', a_i ) 142 zdiag_simin = glob_min( 'icectl', sv_i ) 138 143 zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 139 144 zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) … … 155 160 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rdt_ice 156 161 ! check negative values 157 IF( zdiag_vmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vmin 158 IF( zdiag_amin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_amin 159 IF( zdiag_smin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_smin 160 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 161 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 IF( zdiag_vimin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vimin 163 IF( zdiag_vsmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_s < 0 = ',zdiag_vsmin 164 IF( zdiag_vpmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_ip < 0 = ',zdiag_vpmin 165 IF( zdiag_vlmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_il < 0 = ',zdiag_vlmin 166 IF( zdiag_aimin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_aimin 167 IF( zdiag_simin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_simin 168 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 169 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 170 ! check maximum ice concentration 163 IF( zdiag_a max >MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_a max171 IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 172 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_aimax 165 173 ! check if advection scheme is conservative 166 174 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & … … 193 201 ! water flux 194 202 ! -- mass diag -- ! 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub&196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t )203 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd & 204 & + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 197 205 198 206 ! -- salt diag -- ! … … 200 208 201 209 ! -- heat diag -- ! 202 zdiag_heat 210 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 211 ! equivalent to this: 204 212 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & … … 245 253 IF( icount == 0 ) THEN 246 254 247 pdiag_v = SUM( v_i * rhoi + v_s * rhos , dim=3 )248 pdiag_s = SUM( sv_i * rhoi 255 pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 256 pdiag_s = SUM( sv_i * rhoi , dim=3 ) 249 257 pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 250 258 … … 261 269 262 270 ! -- mass diag -- ! 263 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos , dim=3 ) - pdiag_v ) * r1_rdtice&271 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_rdtice & 264 272 & + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 265 273 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & … … 352 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 353 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 362 ! mean state 363 CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) ! 364 CALL iom_rstput( 0, 0, inum, 'icevol' , SUM(v_i ,dim=3) , ktype = jp_r8 ) ! 365 CALL iom_rstput( 0, 0, inum, 'snwvol' , SUM(v_s ,dim=3) , ktype = jp_r8 ) ! 366 CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) ! 367 CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) ! 354 368 355 369 CALL iom_close( inum ) -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn.F90
r13284 r14026 80 80 ! 81 81 ! controls 82 IF( ln_timing ) CALL timing_start('ice dyn')82 IF( ln_timing ) CALL timing_start('ice_dyn') 83 83 ! 84 84 IF( kt == nit000 .AND. lwp ) THEN … … 175 175 ! 176 176 ! controls 177 IF( ln_timing ) CALL timing_stop ('ice dyn')177 IF( ln_timing ) CALL timing_stop ('ice_dyn') 178 178 ! 179 179 END SUBROUTINE ice_dyn -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90
r13634 r14026 156 156 157 157 ! diagnostics 158 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 158 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 159 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 159 160 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 160 161 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 322 323 ! 323 324 ! --- diagnostics --- ! 324 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 325 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 326 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 325 327 & - zdiag_adv_mass(:,:) ) * z1_dt 326 328 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90
r13617 r14026 183 183 184 184 ! diagnostics 185 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 185 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 186 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 186 187 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 187 188 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 384 385 ! 385 386 ! --- diagnostics --- ! 386 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 387 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 388 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 387 389 & - zdiag_adv_mass(:,:) ) * z1_dt 388 390 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceistate.F90
r13284 r14026 403 403 ! 4) Snow-ice mass (case ice is fully embedded) 404 404 !---------------------------------------------- 405 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s (:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass405 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3 ) ! snow+ice mass 406 406 snwice_mass_b(:,:) = snwice_mass(:,:) 407 407 ! -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceitd.F90
r13617 r14026 29 29 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 30 30 USE prtctl ! Print control 31 USE timing ! Timing 31 32 32 33 IMPLICIT NONE … … 85 86 REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories 86 87 !!------------------------------------------------------------------ 88 IF( ln_timing ) CALL timing_start('iceitd_rem') 87 89 88 90 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' … … 318 320 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 319 321 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 322 IF( ln_timing ) CALL timing_stop ('iceitd_rem') 320 323 ! 321 324 END SUBROUTINE ice_itd_rem … … 481 484 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 482 485 ! 483 ztrans = v_ip_2d(ji,jl1) * zworka(ji) ! Pond volume (also proportional to da/a) 486 !!$ ztrans = v_ip_2d(ji,jl1) * zworka(ji) ! Pond volume (also proportional to da/a) 487 ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume 484 488 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 485 489 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 486 490 ! 487 491 IF ( ln_pnd_lids ) THEN ! Pond lid volume 488 ztrans = v_il_2d(ji,jl1) * zworka(ji) 492 !!$ ztrans = v_il_2d(ji,jl1) * zworka(ji) 493 ztrans = v_il_2d(ji,jl1) * zworkv(ji) 489 494 v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 490 495 v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans … … 592 597 REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred 593 598 !!------------------------------------------------------------------ 599 IF( ln_timing ) CALL timing_start('iceitd_reb') 594 600 ! 595 601 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' … … 623 629 jdonor(ji,jl) = jl 624 630 ! how much of a_i you send in cat sup is somewhat arbitrary 625 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 626 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 627 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 628 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 629 !! zdaice(ji,jl) = a_i_1d(ji) 630 !! zdvice(ji,jl) = v_i_1d(ji) 631 !!clem: these are from UCL and work ok 632 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 633 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 631 ! these are from CICE => transfer everything 632 !!zdaice(ji,jl) = a_i_1d(ji) 633 !!zdvice(ji,jl) = v_i_1d(ji) 634 ! these are from LLN => transfer only half of the category 635 zdaice(ji,jl) = 0.5_wp * a_i_1d(ji) 636 zdvice(ji,jl) = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 634 637 END DO 635 638 ! … … 677 680 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 678 681 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 682 IF( ln_timing ) CALL timing_stop ('iceitd_reb') 679 683 ! 680 684 END SUBROUTINE ice_itd_reb -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icesbc.F90
r13284 r14026 61 61 !!------------------------------------------------------------------- 62 62 ! 63 IF( ln_timing ) CALL timing_start('ice _sbc')63 IF( ln_timing ) CALL timing_start('icesbc') 64 64 ! 65 65 IF( kt == nit000 .AND. lwp ) THEN … … 86 86 ENDIF 87 87 ! 88 IF( ln_timing ) CALL timing_stop('ice _sbc')88 IF( ln_timing ) CALL timing_stop('icesbc') 89 89 ! 90 90 END SUBROUTINE ice_sbc_tau … … 119 119 !!-------------------------------------------------------------------- 120 120 ! 121 IF( ln_timing ) CALL timing_start('ice _sbc_flx')121 IF( ln_timing ) CALL timing_start('icesbc') 122 122 123 123 IF( kt == nit000 .AND. lwp ) THEN … … 172 172 ENDIF 173 173 ! 174 IF( ln_timing ) CALL timing_stop('ice _sbc_flx')174 IF( ln_timing ) CALL timing_stop('icesbc') 175 175 ! 176 176 END SUBROUTINE ice_sbc_flx -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90
r13720 r14026 120 120 !!---------------------------------------------------------------------- 121 121 ! 122 IF( ln_timing ) CALL timing_start('ice _stp')122 IF( ln_timing ) CALL timing_start('icestp') 123 123 ! 124 124 ! !-----------------------! … … 214 214 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 215 215 ! 216 IF( ln_timing ) CALL timing_stop('ice _stp')216 IF( ln_timing ) CALL timing_stop('icestp') 217 217 ! 218 218 END SUBROUTINE ice_stp … … 365 365 v_i_b (:,:,:) = v_i (:,:,:) ! ice volume 366 366 v_s_b (:,:,:) = v_s (:,:,:) ! snow volume 367 v_ip_b(:,:,:) = v_ip(:,:,:) ! pond volume 368 v_il_b(:,:,:) = v_il(:,:,:) ! pond lid volume 367 369 sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content 368 370 e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy … … 425 427 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 426 428 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 427 diag_aice(ji,jj) = 0._wp 429 diag_aice(ji,jj) = 0._wp ; diag_vpnd(ji,jj) = 0._wp 428 430 429 431 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) … … 482 484 diag_vsnw(:,:) = diag_vsnw(:,:) & 483 485 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_rdtice * rhos 486 diag_vpnd(:,:) = diag_vpnd(:,:) & 487 & + SUM( v_ip + v_il - v_ip_b - v_il_b , dim=3 ) * r1_rdtice * rhow 484 488 ! 485 489 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90
r13642 r14026 171 171 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 172 172 173 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 174 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 175 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 176 zqfr = 0._wp 177 zqfr_pos = 0._wp 178 qsb_ice_bot(ji,jj) = 0._wp 179 ENDIF 180 ! 173 181 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 174 182 ! qlead is the energy received from the atm. in the leads. … … 247 255 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 248 256 CALL ice_thd_dh ! Ice-Snow thickness 249 CALL ice_thd_pnd ! Melt ponds formation250 257 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 251 258 ENDIF … … 258 265 ! 259 266 IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- ! 267 ! 268 IF( ln_pnd .AND. ln_icedH ) & 269 & CALL ice_thd_pnd ! --- Melt ponds formation --- ! 260 270 ! 261 271 CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_dh.F90
r13642 r14026 55 55 !! - Snow ice formation 56 56 !! 57 !! ** Note : h=max(0,h+dh) are often used to ensure positivity of h. 58 !! very small negative values can occur otherwise (e.g. -1.e-20) 59 !! 57 60 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 58 61 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 … … 79 82 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 80 83 81 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)82 84 REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2) 83 85 REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2) … … 85 87 REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 86 88 REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2) 87 88 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 89 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 90 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 91 92 REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness 93 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 94 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 95 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 96 89 REAL(wp), DIMENSION(jpij) :: zdeltah 97 90 REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing 91 92 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanishing by melting 93 REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: zh_i ! ice layer thickness (m) 94 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: zh_s ! snw layer thickness (m) 95 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3) 98 96 99 97 REAL(wp) :: zswitch_sal … … 108 106 END SELECT 109 107 110 ! initialize layer thicknesses and enthalpies 108 ! initialize ice layer thicknesses and enthalpies 109 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 111 110 h_i_old (1:npti,0:nlay_i+1) = 0._wp 112 eh_i_old(1:npti,0:nlay_i+1) = 0._wp111 zh_i (1:npti,0:nlay_i+1) = 0._wp 113 112 DO jk = 1, nlay_i 114 113 DO ji = 1, npti 114 eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 115 115 h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 116 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 116 zh_i (ji,jk) = h_i_1d(ji) * r1_nlay_i 117 END DO 118 END DO 119 ! 120 ! initialize snw layer thicknesses and enthalpies 121 zh_s(1:npti,0) = 0._wp 122 ze_s(1:npti,0) = 0._wp 123 DO jk = 1, nlay_s 124 DO ji = 1, npti 125 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 126 ze_s(ji,jk) = e_s_1d(ji,jk) 117 127 END DO 118 128 END DO … … 141 151 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 142 152 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 143 END DO144 145 ! Ice and snow layer thicknesses146 !-------------------------------147 DO jk = 1, nlay_i148 DO ji = 1, npti149 zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i150 END DO151 END DO152 153 DO jk = 1, nlay_s154 DO ji = 1, npti155 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s156 END DO157 153 END DO 158 154 … … 167 163 DO ji = 1, npti 168 164 IF( t_s_1d(ji,jk) > rt0 ) THEN 169 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0170 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos 165 hfx_res_1d (ji) = hfx_res_1d (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0 166 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux 171 167 ! updates 172 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)173 h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)168 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) 169 h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(ji,jk) ) 174 170 zh_s (ji,jk) = 0._wp 175 e_s_1d (ji,jk) = 0._wp 176 t_s_1d (ji,jk) = rt0 171 ze_s (ji,jk) = 0._wp 177 172 END IF 178 173 END DO … … 181 176 ! Snow precipitation 182 177 !------------------- 183 CALL ice_var_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 184 185 zdeltah(1:npti,:) = 0._wp 178 CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 179 186 180 DO ji = 1, npti 187 181 IF( sprecip_1d(ji) > 0._wp ) THEN 182 zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip 183 ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3) 188 184 ! 189 ! --- precipitation --- 190 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji) ! thickness change 191 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 185 hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_rdtice ! heat flux from snow precip (>0, W.m-2) 186 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * zh_s(ji,0) * a_i_1d(ji) * r1_rdtice ! mass flux, <0 192 187 ! 193 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat flux from snow precip (>0, W.m-2)194 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! mass flux, <0195 196 ! --- melt of falling snow ---197 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )198 zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change199 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting200 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat used to melt snow (W.m-2, >0)201 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0202 203 ! updates available heat + precipitations after melting204 dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)205 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )206 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)207 208 188 ! update thickness 209 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 210 ! 211 ELSE 212 ! 213 zdh_s_pre(ji) = 0._wp 214 zqprec (ji) = 0._wp 215 ! 189 h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 216 190 ENDIF 217 END DO218 219 ! recalculate snow layers220 DO jk = 1, nlay_s221 DO ji = 1, npti222 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s223 END DO224 191 END DO 225 192 226 193 ! Snow melting 227 194 ! ------------ 228 ! If heat still available (zq_top > 0), then melt more snow 229 zdeltah(1:npti,:) = 0._wp 230 zdh_s_mel(1:npti) = 0._wp 231 DO jk = 1, nlay_s 195 ! If heat still available (zq_top > 0) 196 ! then all snw precip has been melted and we need to melt more snow 197 DO jk = 0, nlay_s 232 198 DO ji = 1, npti 233 199 IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 234 200 ! 235 rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 236 zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 237 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 238 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 239 240 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) 241 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) 201 rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) ) 202 zdum = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 ) ! thickness change 203 zdum = MAX( zdum , - zh_s(ji,jk) ) ! bound melting 204 205 hfx_snw_1d (ji) = hfx_snw_1d (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice ! heat used to melt snow(W.m-2, >0) 206 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_rdtice ! snow melting only = water into the ocean 242 207 243 208 ! updates available heat + thickness 244 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 245 zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 246 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 247 zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 209 dh_s_mlt(ji) = dh_s_mlt(ji) + zdum 210 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(ji,jk) ) 211 h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum ) 212 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 213 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 248 214 ! 249 215 ENDIF … … 255 221 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 256 222 ! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 257 zdeltah(1:npti,:) = 0._wp 223 zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 224 zevap_rema(1:npti) = 0._wp 258 225 DO ji = 1, npti 259 226 IF( evap_ice_1d(ji) > 0._wp ) THEN 227 zdeltah (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rdt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, < 0 228 zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rdt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on) 229 ENDIF 230 END DO 231 232 DO jk = 0, nlay_s 233 DO ji = 1, npti 234 zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 260 235 ! 261 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice ) 262 zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos ! remaining evap in kg.m-2 (used for ice melting later on) 263 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 264 265 hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 266 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 267 & * a_i_1d(ji) * r1_rdtice 268 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! Mass flux by sublimation 269 270 ! new snow thickness 271 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 272 ! update precipitations after sublimation and correct sublimation 273 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 274 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 275 ! 276 ELSE 277 ! 278 zdh_s_sub (ji) = 0._wp 279 zevap_rema(ji) = 0._wp 280 ! 281 ENDIF 282 END DO 283 284 ! --- Update snow diags --- ! 285 DO ji = 1, npti 286 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 287 END DO 288 289 ! Update temperature, energy 290 !--------------------------- 291 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 292 DO jk = 1, nlay_s 293 DO ji = 1,npti 294 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 295 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 296 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 297 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 298 END DO 299 END DO 300 236 hfx_sub_1d (ji) = hfx_sub_1d (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice ! Heat flux of snw that sublimates [W.m-2], < 0 237 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux by sublimation 238 239 ! update thickness 240 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum ) 241 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 242 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 243 244 ! update sublimation left 245 zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) 246 END DO 247 END DO 248 249 ! 301 250 ! ! ============ ! 302 251 ! ! Ice ! … … 305 254 ! Surface ice melting 306 255 !-------------------- 307 zdeltah(1:npti,:) = 0._wp ! important308 256 DO jk = 1, nlay_i 309 257 DO ji = 1, npti … … 313 261 314 262 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 315 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 316 ! set up at 0 since no energy is needed to melt water...(it is already melted) 317 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 318 ! this should normally not happen, but sometimes, heat diffusion leads to this 319 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 320 321 dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting 322 323 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 324 325 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 326 ! ice enthalpy zEi is "sent" to the ocean 327 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 328 ! using s_i_1d and not sz_i_1d(jk) is ok 329 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 330 263 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 264 ! set up at 0 since no energy is needed to melt water...(it is already melted) 265 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 266 ! this should normally not happen, but sometimes, heat diffusion leads to this 267 zfmdt = - zdum * rhoi ! Recompute mass flux [kg/m2, >0] 268 ! 269 dh_i_itm(ji) = dh_i_itm(ji) + zdum ! Cumulate internal melting 270 ! 271 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 272 ! ice enthalpy zEi is "sent" to the ocean 273 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux 274 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux 275 ! using s_i_1d and not sz_i_1d(jk) is ok 331 276 ELSE !-- Surface melting 332 277 … … 337 282 zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 338 283 339 zd eltah(ji,jk)= - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]340 341 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]342 343 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk)* rhoi * zdE ) ! update available heat344 345 dh_i_sum(ji) = dh_i_sum(ji) + zd eltah(ji,jk)! Cumulate surface melt346 347 zfmdt = - rhoi * zd eltah(ji,jk)! Recompute mass flux [kg/m2, >0]284 zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0] 285 286 zdum = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 287 288 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat 289 290 dh_i_sum(ji) = dh_i_sum(ji) + zdum ! Cumulate surface melt 291 292 zfmdt = - rhoi * zdum ! Recompute mass flux [kg/m2, >0] 348 293 349 294 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 350 295 351 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux >0 352 ! using s_i_1d and not sz_i_1d(jk) is ok) 353 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux [W.m-2], < 0 354 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 355 ! 356 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 357 296 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0 297 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux used in this process [W.m-2], > 0 298 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux 299 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux >0 300 ! using s_i_1d and not sz_i_1d(jk) is ok) 358 301 END IF 359 302 ! update thickness 303 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 304 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 305 ! 306 ! update heat content (J.m-2) and layer thickness 307 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 308 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 309 ! 310 ! 360 311 ! Ice sublimation 361 312 ! --------------- 362 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 363 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 364 dh_i_sub(ji) = dh_i_sub(ji) + zdum 365 366 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 367 ! clem: flux is sent to the ocean for simplicity 368 ! but salt should remain in the ice except 369 ! if all ice is melted. => must be corrected 370 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0 371 372 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice ! Mass flux > 0 373 374 ! update remaining mass flux 375 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 376 313 zdum = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) 314 ! 315 hfx_sub_1d(ji) = hfx_sub_1d(ji) + e_i_1d(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0 316 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux > 0 317 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux >0 318 ! clem: flux is sent to the ocean for simplicity 319 ! but salt should remain in the ice except 320 ! if all ice is melted. => must be corrected 321 ! update remaining mass flux and thickness 322 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 323 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 324 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 325 dh_i_sub(ji) = dh_i_sub(ji) + zdum 326 327 ! update heat content (J.m-2) and layer thickness 328 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 329 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 330 377 331 ! record which layers have disappeared (for bottom melting) 378 332 ! => icount=0 : no layer has vanished 379 333 ! => icount=5 : 5 layers have vanished 380 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk)) ) )334 rswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) ) 381 335 icount(ji,jk) = NINT( rswitch ) 382 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )383 336 384 ! update heat content (J.m-2) and layer thickness 385 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 386 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 387 END DO 388 END DO 389 390 ! update ice thickness 391 DO ji = 1, npti 392 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 393 END DO 394 337 END DO 338 END DO 339 395 340 ! remaining "potential" evap is sent to ocean 396 341 DO ji = 1, npti … … 430 375 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 431 376 432 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity433 434 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)435 436 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)437 438 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)439 & - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp* ztmelts440 441 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)442 443 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)444 445 dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )377 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity 378 379 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C) 380 381 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 382 383 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 384 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 385 386 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 387 388 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 389 390 dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 446 391 447 392 END DO 448 393 ! Contribution to Energy and Salt Fluxes 449 zfmdt = - rhoi * dh_i_bog(ji)! Mass flux x time step (kg/m2, < 0)394 zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0) 450 395 451 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], >0 452 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], <0 453 454 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice ! Salt flux, <0 455 456 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice ! Mass flux, <0 396 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux to the ocean [W.m-2], >0 397 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux used in this process [W.m-2], <0 398 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * a_i_1d(ji) * r1_rdtice ! Mass flux, <0 399 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux, <0 400 401 ! update thickness 402 zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji) 403 h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji) 457 404 458 405 ! update heat content (J.m-2) and layer thickness … … 466 413 ! Ice Basal melt 467 414 !--------------- 468 zdeltah(1:npti,:) = 0._wp ! important469 415 DO jk = nlay_i, 1, -1 470 416 DO ji = 1, npti … … 475 421 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 476 422 477 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 478 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 479 ! set up at 0 since no energy is needed to melt water...(it is already melted) 480 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 481 ! this should normally not happen, but sometimes, heat diffusion leads to this 482 483 dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk) 484 485 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 486 487 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 488 ! ice enthalpy zEi is "sent" to the ocean 489 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 490 ! using s_i_1d and not sz_i_1d(jk) is ok 491 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 492 493 ! update heat content (J.m-2) and layer thickness 494 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 495 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 496 423 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 424 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 425 ! set up at 0 since no energy is needed to melt water...(it is already melted) 426 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 427 ! this should normally not happen, but sometimes, heat diffusion leads to this 428 dh_i_itm (ji) = dh_i_itm(ji) + zdum 429 ! 430 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 431 ! 432 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 433 ! ice enthalpy zEi is "sent" to the ocean 434 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux 435 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux 436 ! using s_i_1d and not sz_i_1d(jk) is ok 497 437 ELSE !-- Basal melting 498 438 499 zEi = - e_i_1d(ji,jk) * r1_rhoi! Specific enthalpy of melting ice (J/kg, <0)500 zEw = rcp * ztmelts! Specific enthalpy of meltwater (J/kg, <0)501 zdE = zEi - zEw! Specific enthalpy difference (J/kg, <0)502 503 zfmdt = - zq_bot(ji) / zdE! Mass flux x time step (kg/m2, >0)504 505 zd eltah(ji,jk) = - zfmdt * r1_rhoi! Gross thickness change506 507 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change439 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 440 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 441 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 442 443 zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0) 444 445 zdum = - zfmdt * r1_rhoi ! Gross thickness change 446 447 zdum = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) ) ! bound thickness change 508 448 509 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 510 511 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt 512 513 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 514 515 zQm = zfmdt * zEw ! Heat exchanged with ocean 516 517 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 518 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat used in this process [W.m-2], >0 519 520 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 521 ! using s_i_1d and not sz_i_1d(jk) is ok 522 523 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 524 525 ! update heat content (J.m-2) and layer thickness 526 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 527 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 449 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 450 451 dh_i_bom(ji) = dh_i_bom(ji) + zdum ! Update basal melt 452 453 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 454 455 zQm = zfmdt * zEw ! Heat exchanged with ocean 456 457 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 458 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat used in this process [W.m-2], >0 459 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_rdtice ! Mass flux 460 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux 461 ! using s_i_1d and not sz_i_1d(jk) is ok 528 462 ENDIF 529 463 ! update thickness 464 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 465 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 466 ! 467 ! update heat content (J.m-2) and layer thickness 468 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 469 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 530 470 ENDIF 531 471 END DO 532 472 END DO 533 473 534 ! Update temperature, energy 535 ! -------------------------- 536 DO ji = 1, npti 537 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 538 END DO 539 540 ! If heat still available then melt more snow 541 !------------------------------------------- 542 zdeltah(1:npti,:) = 0._wp ! important 543 DO ji = 1, npti 544 zq_rema (ji) = zq_top(ji) + zq_bot(ji) 545 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow 546 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 547 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 548 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 549 dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1) 550 h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1) 551 552 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 553 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! Heat used to melt snow, W.m-2 (>0) 554 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! Mass flux 555 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 556 ! 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 !!!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 559 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 561 END DO 562 563 ! 474 ! Remove snow if ice has melted entirely 475 ! -------------------------------------- 476 DO jk = 0, nlay_s 477 DO ji = 1,npti 478 IF( h_i_1d(ji) == 0._wp ) THEN 479 ! mass & energy loss to the ocean 480 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0 481 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux 482 483 ! update thickness and energy 484 h_s_1d(ji) = 0._wp 485 ze_s (ji,jk) = 0._wp 486 zh_s (ji,jk) = 0._wp 487 ENDIF 488 END DO 489 END DO 490 491 ! Snow load on ice 492 ! ----------------- 493 ! When snow load exceeds Archimede's limit and sst is positive, 494 ! snow-ice formation (next bloc) can lead to negative ice enthalpy. 495 ! Therefore we consider here that this excess of snow falls into the ocean 496 zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rau0) * r1_rhos 497 DO jk = 0, nlay_s 498 DO ji = 1, npti 499 IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN 500 ! snow layer thickness that falls into the ocean 501 zdum = MIN( zdeltah(ji) , zh_s(ji,jk) ) 502 ! mass & energy loss to the ocean 503 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0 504 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_rdtice ! mass flux 505 ! update thickness and energy 506 h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum ) 507 zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) 508 ! update snow thickness that still has to fall 509 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 510 ENDIF 511 END DO 512 END DO 513 564 514 ! Snow-Ice formation 565 515 ! ------------------ 566 ! When snow load exce sses Archimede's limit, snow-ice interface goes down under sea-level,567 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice516 ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level, 517 ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice) 568 518 z1_rho = 1._wp / ( rhos+rau0-rhoi ) 519 zdeltah(1:npti) = 0._wp 569 520 DO ji = 1, npti 570 521 ! 571 dh_snowice(ji) = MAX( 522 dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho ) 572 523 573 524 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 579 530 zQm = zfmdt * zEw 580 531 581 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux 582 583 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux 532 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux 533 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_rdtice ! Salt flux 584 534 585 535 ! Case constant salinity in time: virtual salt flux to keep salinity constant 586 536 IF( nn_icesal /= 2 ) THEN 587 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_rdtice &! put back sss_m into the ocean588 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi* r1_rdtice ! and get rn_icesal from the ocean537 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_rdtice & ! put back sss_m into the ocean 538 & - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_rdtice ! and get rn_icesal from the ocean 589 539 ENDIF 590 540 591 541 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 592 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice 593 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice 542 wfx_sni_1d (ji) = wfx_sni_1d (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_rdtice 543 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_rdtice 544 545 ! update thickness 546 zh_i(ji,0) = zh_i(ji,0) + dh_snowice(ji) 547 zdeltah(ji) = dh_snowice(ji) 594 548 595 549 ! update heat content (J.m-2) and layer thickness 596 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw597 550 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 598 599 END DO 600 601 ! 602 ! Update temperature, energy 603 ! -------------------------- 604 DO ji = 1, npti 605 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 606 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 607 END DO 608 551 eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw ! 1st part (sea water enthalpy) 552 553 END DO 554 ! 555 DO jk = nlay_s, 0, -1 ! flooding of snow starts from the base 556 DO ji = 1, npti 557 zdum = MIN( zdeltah(ji), zh_s(ji,jk) ) ! amount of snw that floods, > 0 558 zh_s(ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) ! remove some snow thickness 559 eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy) 560 ! update dh_snowice 561 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 562 END DO 563 END DO 564 ! 565 ! 566 !!$ ! --- Update snow diags --- ! 567 !!$ !!clem: this is wrong. dh_s_tot is not used anyway 568 !!$ DO ji = 1, npti 569 !!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji) 570 !!$ END DO 571 ! 572 ! 573 ! Remapping of snw enthalpy on a regular grid 574 !-------------------------------------------- 575 CALL snw_ent( zh_s, ze_s, e_s_1d ) 576 577 ! recalculate t_s_1d from e_s_1d 609 578 DO jk = 1, nlay_s 610 579 DO ji = 1,npti 611 ! where there is no ice or no snow 612 rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) ) 613 ! mass & energy loss to the ocean 614 hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 615 & ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice ) ! heat flux to the ocean [W.m-2], < 0 616 wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 617 & ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice ) ! mass flux 618 ! update energy (mass is updated in the next loop) 619 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 620 ! recalculate t_s_1d from e_s_1d 621 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 622 END DO 623 END DO 580 IF( h_s_1d(ji) > 0._wp ) THEN 581 t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 582 ELSE 583 t_s_1d(ji,jk) = rt0 584 ENDIF 585 END DO 586 END DO 587 588 ! Note: remapping of ice enthalpy is done in icethd.F90 624 589 625 590 ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 626 591 WHERE( h_i_1d(1:npti) == 0._wp ) 627 a_i_1d(1:npti) = 0._wp 628 h_s_1d(1:npti) = 0._wp 592 a_i_1d (1:npti) = 0._wp 593 h_s_1d (1:npti) = 0._wp 594 t_su_1d(1:npti) = rt0 629 595 END WHERE 630 !596 631 597 END SUBROUTINE ice_thd_dh 632 598 599 SUBROUTINE snw_ent( ph_old, pe_old, pe_new ) 600 !!------------------------------------------------------------------- 601 !! *** ROUTINE snw_ent *** 602 !! 603 !! ** Purpose : 604 !! This routine computes new vertical grids in the snow, 605 !! and consistently redistributes temperatures. 606 !! Redistribution is made so as to ensure to energy conservation 607 !! 608 !! 609 !! ** Method : linear conservative remapping 610 !! 611 !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 612 !! 2) linear remapping on the new layers 613 !! 614 !! ------------ cum0(0) ------------- cum1(0) 615 !! NEW ------------- 616 !! ------------ cum0(1) ==> ------------- 617 !! ... ------------- 618 !! ------------ ------------- 619 !! ------------ cum0(nlay_s+1) ------------- cum1(nlay_s) 620 !! 621 !! 622 !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 623 !!------------------------------------------------------------------- 624 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: ph_old ! old thicknesses (m) 625 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: pe_old ! old enthlapies (J.m-3) 626 REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped) 627 ! 628 INTEGER :: ji ! dummy loop indices 629 INTEGER :: jk0, jk1 ! old/new layer indices 630 ! 631 REAL(wp), DIMENSION(jpij,0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces 632 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces 633 REAL(wp), DIMENSION(jpij) :: zhnew ! new layers thicknesses 634 !!------------------------------------------------------------------- 635 636 !-------------------------------------------------------------------------- 637 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 638 !-------------------------------------------------------------------------- 639 zeh_cum0(1:npti,0) = 0._wp 640 zh_cum0 (1:npti,0) = 0._wp 641 DO jk0 = 1, nlay_s+1 642 DO ji = 1, npti 643 zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1) 644 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1) 645 END DO 646 END DO 647 648 !------------------------------------ 649 ! 2) Interpolation on the new layers 650 !------------------------------------ 651 ! new layer thickesses 652 DO ji = 1, npti 653 zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s 654 END DO 655 656 ! new layers interfaces 657 zh_cum1(1:npti,0) = 0._wp 658 DO jk1 = 1, nlay_s 659 DO ji = 1, npti 660 zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 661 END DO 662 END DO 663 664 zeh_cum1(1:npti,0:nlay_s) = 0._wp 665 ! new cumulative q*h => linear interpolation 666 DO jk0 = 1, nlay_s+1 667 DO jk1 = 1, nlay_s-1 668 DO ji = 1, npti 669 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 670 zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & 671 & zeh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) & 672 & / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 673 ENDIF 674 END DO 675 END DO 676 END DO 677 ! to ensure that total heat content is strictly conserved, set: 678 zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1) 679 680 ! new enthalpies 681 DO jk1 = 1, nlay_s 682 DO ji = 1, npti 683 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 684 pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 685 END DO 686 END DO 687 688 END SUBROUTINE snw_ent 689 690 633 691 #else 634 692 !!---------------------------------------------------------------------- -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_pnd.F90
r13284 r14026 83 83 !!------------------------------------------------------------------- 84 84 INTEGER :: ji ! loop indices 85 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 85 86 !!------------------------------------------------------------------- 86 87 DO ji = 1, npti 87 88 ! 88 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 89 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 90 ! 91 IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 89 92 h_ip_1d(ji) = rn_hpnd 90 93 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) … … 96 99 ENDIF 97 100 ! 101 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 102 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice 103 ! 98 104 END DO 99 105 ! … … 134 140 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 135 141 !! 136 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi 137 !! perm = permability of sea-ice 142 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi * flush --- from Flocco et al 2007 --- 143 !! perm = permability of sea-ice + correction from Hunke et al 2012 (flush) 138 144 !! visc = water viscosity 139 145 !! Hp = height of top of the pond above sea-level 140 146 !! Hi = ice thickness thru which there is flushing 147 !! flush= correction otherwise flushing is excessive 141 148 !! 142 149 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness … … 147 154 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 148 155 !! 149 !! ** Note : mostly stolen from CICE 156 !! ** Note : mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 150 157 !! 151 158 !! ** References : Flocco and Feltham (JGR, 2007) 152 159 !! Flocco et al (JGR, 2010) 153 160 !! Holland et al (J. Clim, 2012) 161 !! Hunke et al (OM 2012) 154 162 !!------------------------------------------------------------------- 155 163 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array … … 158 166 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 159 167 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 160 !! 161 REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding 168 REAL(wp), PARAMETER :: zflush = 1.e-02_wp ! tuning param to reduce flushing (from Hunke et al 2012) 169 !! 170 REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding 162 171 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 172 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 163 173 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 164 174 REAL(wp) :: zv_ip_max ! max pond volume allowed 165 175 REAL(wp) :: zdT ! zTp-t_su 166 REAL(wp) :: zsbr 176 REAL(wp) :: zsbr, ztmelts ! Brine salinity 167 177 REAL(wp) :: zperm ! permeability of sea ice 168 178 REAL(wp) :: zfac, zdum ! temporary arrays … … 176 186 177 187 DO ji = 1, npti 188 ! 189 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 178 190 ! !----------------------------------------------------! 179 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN! Case ice thickness < rn_himin or tiny ice fraction !191 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 180 192 ! !----------------------------------------------------! 181 193 !--- Remove ponds on thin ice or tiny ice fractions 182 a_ip_1d(ji) 183 h_ip_1d(ji) 184 h_il_1d(ji) 194 a_ip_1d(ji) = 0._wp 195 h_ip_1d(ji) = 0._wp 196 h_il_1d(ji) = 0._wp 185 197 ! !--------------------------------! 186 198 ELSE ! Case ice thickness >= rn_himin ! … … 193 205 !------------------! 194 206 ! 195 !--- available meltwater for melt ponding ---!196 zd um = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)197 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed198 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum) ! max for roundoff errors?207 !--- available meltwater for melt ponding (zdv_avail) ---! 208 zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 209 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed 210 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 199 211 ! 200 212 !--- overflow ---! 201 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 202 ! a_ip_max = zfr_mlt * a_i 203 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 213 ! 214 ! area driven overflow 215 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 216 ! a_ip_max = zfr_mlt * a_i 217 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 204 218 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 205 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 206 207 ! If pond depth exceeds half the ice thickness then reduce the pond volume 208 ! h_ip_max = 0.5 * h_i 209 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 219 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 220 221 ! depth driven overflow 222 ! If pond depth exceeds half the ice thickness then reduce the pond volume 223 ! h_ip_max = 0.5 * h_i 224 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 210 225 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 211 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )226 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 212 227 213 228 !--- Pond growing ---! … … 217 232 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 218 233 ! 219 !--- mass flux ---!220 IF( zdv_mlt > 0._wp ) THEN221 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1]222 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac223 !224 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux225 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)226 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum)227 ENDIF228 229 234 !-------------------! 230 235 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) … … 237 242 ! 238 243 !--- Lid growing and subsequent pond shrinking ---! 239 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0244 zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 240 245 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 241 246 242 247 ! Lid growing 243 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) +zdv_frz )248 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 244 249 245 250 ! Pond shrinking 246 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) -zdv_frz )251 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 247 252 248 253 ELSE 254 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 249 255 ! Pond shrinking 250 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)256 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 251 257 ENDIF 252 258 ! … … 254 260 ! v_ip = h_ip * a_ip 255 261 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 256 a_ip_1d(ji) 257 h_ip_1d(ji) 258 259 !--------------- !260 ! Pond flushing!261 !--------------- !262 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 263 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 264 265 !------------------------------------------------! 266 ! Pond drainage through brine network (flushing) ! 267 !------------------------------------------------! 262 268 ! height of top of the pond above sea-level 263 269 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 … … 265 271 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 266 272 DO jk = 1, nlay_i 267 zsbr = - 1.2_wp & 268 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 269 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 270 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 271 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 273 ! MV Assur is inconsistent with SI3 274 !!zsbr = - 1.2_wp & 275 !! & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 276 !! & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 277 !! & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 278 !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 279 ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 280 ztmelts = -rTmlt * sz_i_1d(ji,jk) 281 ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 272 282 END DO 273 283 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 274 284 275 285 ! Do the drainage using Darcy's law 276 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 277 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 286 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * zflush ! zflush comes from Hunke et al. (2012) 287 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 278 288 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 279 289 280 290 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 281 a_ip_1d(ji) 282 h_ip_1d(ji) 291 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 292 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 283 293 284 294 !--- Corrections and lid thickness ---! 285 295 IF( ln_pnd_lids ) THEN 286 296 !--- retrieve lid thickness from volume ---! 287 IF( a_ip_1d(ji) > epsi10) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)288 ELSE ; h_il_1d(ji) = 0._wp297 IF( a_ip_1d(ji) > 0.01_wp ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 298 ELSE ; h_il_1d(ji) = 0._wp 289 299 ENDIF 290 300 !--- remove ponds if lids are much larger than ponds ---! 291 301 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 292 a_ip_1d(ji) 293 h_ip_1d(ji) 294 h_il_1d(ji) 302 a_ip_1d(ji) = 0._wp 303 h_ip_1d(ji) = 0._wp 304 h_il_1d(ji) = 0._wp 295 305 ENDIF 296 306 ENDIF 307 308 !!$ ! diagnostics: dvpnd = mlt+rnf+frz+drn 309 !!$ diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow * zdv_avail * r1_rdtice ! > 0, surface melt input 310 !!$ diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_rdtice ! < 0, runoff 311 !!$ diag_dvpn_frz_1d(ji) = diag_dvpn_frz_1d(ji) + rhow * zdv_frz * r1_rdtice ! < 0, shrinking 312 !!$ diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow * zdv_flush * r1_rdtice ! < 0, drainage 297 313 ! 298 314 ENDIF 299 315 ! 316 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 317 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice 318 ! 300 319 END DO 301 320 ! -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_zdf_bl99.F90
r13284 r14026 109 109 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 110 110 ! 111 REAL(wp), DIMENSION(jpij ) 112 REAL(wp), DIMENSION(jpij,nlay_i) 113 REAL(wp), DIMENSION(jpij,nlay_s) 114 REAL(wp), DIMENSION(jpij,nlay_i) 115 REAL(wp), DIMENSION(jpij,nlay_s) 116 REAL(wp), DIMENSION(jpij,0:nlay_i) 117 REAL(wp), DIMENSION(jpij,0:nlay_i) 118 REAL(wp), DIMENSION(jpij,0:nlay_i) 119 REAL(wp), DIMENSION(jpij,0:nlay_i) 120 REAL(wp), DIMENSION(jpij,0:nlay_i) 121 REAL(wp), DIMENSION(jpij,0:nlay_i) 122 REAL(wp), DIMENSION(jpij,0:nlay_s) 123 REAL(wp), DIMENSION(jpij,0:nlay_s) 124 REAL(wp), DIMENSION(jpij,0:nlay_s) 125 REAL(wp), DIMENSION(jpij,0:nlay_s) 126 REAL(wp), DIMENSION(jpij) 127 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindterm ! 'Ind'ependent term128 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term129 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term130 REAL(wp), DIMENSION(jpij ,nlay_i+3,3) :: ztrid ! Tridiagonal system terms131 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat132 REAL(wp), DIMENSION(jpij ) :: zghe ! G(he), th. conduct enhancement factor, mono-cat133 REAL(wp), DIMENSION(jpij ) :: za_s_fra ! ice fraction covered by snow134 REAL(wp), DIMENSION(jpij ) :: isnow ! snow presence (1) or not (0)135 REAL(wp), DIMENSION(jpij ) :: isnow_comb ! snow presence for met-office111 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 112 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 113 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow 114 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 122 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 123 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 124 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 126 REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity 127 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 128 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 129 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 132 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term 133 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term 134 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zdiagbis ! Temporary 'dia'gonal term 135 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) :: ztrid ! Tridiagonal system terms 136 136 ! 137 137 ! Mono-category … … 533 533 ! Solve the tridiagonal system with Gauss elimination method. 534 534 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 535 jm_maxt = 0 536 jm_mint = nlay_i+5 537 DO ji = 1, npti 538 jm_mint = MIN(jm_min(ji),jm_mint) 539 jm_maxt = MAX(jm_max(ji),jm_maxt) 540 END DO 541 542 DO jk = jm_mint+1, jm_maxt 543 DO ji = 1, npti 544 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 535 !!$ jm_maxt = 0 536 !!$ jm_mint = nlay_i+5 537 !!$ DO ji = 1, npti 538 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 539 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 540 !!$ END DO 541 !!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 542 !!$ 543 !!$ DO jk = jm_mint+1, jm_maxt 544 !!$ DO ji = 1, npti 545 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 546 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 547 !!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 548 !!$ END DO 549 !!$ END DO 550 ! clem: maybe one should find a way to reverse this loop for mpi performance 551 DO ji = 1, npti 552 jm_mint = jm_min(ji) 553 jm_maxt = jm_max(ji) 554 DO jm = jm_mint+1, jm_maxt 545 555 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 546 556 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) … … 564 574 END DO 565 575 576 ! snow temperatures 566 577 DO ji = 1, npti 567 578 ! Variables used after iterations 568 579 ! Value must be frozen after convergence for MPP independance reason 569 IF ( .NOT. l_T_converged(ji) ) THEN 570 ! snow temperatures 571 IF( h_s_1d(ji) > 0._wp ) THEN 572 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 573 ENDIF 574 ! surface temperature 580 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 581 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 582 END DO 583 !!clem SNWLAY 584 DO jm = nlay_s, 2, -1 585 DO ji = 1, npti 586 jk = jm - 1 587 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 588 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 589 END DO 590 END DO 591 592 ! surface temperature 593 DO ji = 1, npti 594 IF( .NOT. l_T_converged(ji) ) THEN 575 595 ztsub(ji) = t_su_1d(ji) 576 596 IF( t_su_1d(ji) < rt0 ) THEN 577 t_su_1d(ji) = ( 578 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) *t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))597 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 598 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 579 599 ENDIF 580 600 ENDIF 581 601 END DO 582 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)583 602 ! 584 603 !-------------------------------------------------------------- … … 727 746 ! Solve the tridiagonal system with Gauss elimination method. 728 747 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 729 jm_maxt = 0 730 jm_mint = nlay_i+5 731 DO ji = 1, npti 732 jm_mint = MIN(jm_min(ji),jm_mint) 733 jm_maxt = MAX(jm_max(ji),jm_maxt) 734 END DO 735 736 DO jk = jm_mint+1, jm_maxt 737 DO ji = 1, npti 738 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 748 !!$ jm_maxt = 0 749 !!$ jm_mint = nlay_i+5 750 !!$ DO ji = 1, npti 751 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 752 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 753 !!$ END DO 754 !!$ 755 !!$ DO jk = jm_mint+1, jm_maxt 756 !!$ DO ji = 1, npti 757 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 758 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 759 !!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 760 !!$ END DO 761 !!$ END DO 762 ! clem: maybe one should find a way to reverse this loop for mpi performance 763 DO ji = 1, npti 764 jm_mint = jm_min(ji) 765 jm_maxt = jm_max(ji) 766 DO jm = jm_mint+1, jm_maxt 739 767 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 740 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1)/ zdiagbis(ji,jm-1)768 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 741 769 END DO 742 770 END DO 743 771 744 772 ! ice temperatures 745 773 DO ji = 1, npti … … 761 789 ! snow temperatures 762 790 DO ji = 1, npti 763 ! Variable used after iterations791 ! Variables used after iterations 764 792 ! Value must be frozen after convergence for MPP independance reason 765 IF ( .NOT. l_T_converged(ji) ) THEN 766 IF( h_s_1d(ji) > 0._wp ) THEN 767 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 768 ENDIF 769 ENDIF 770 END DO 771 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 793 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 794 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 795 END DO 796 !!clem SNWLAY 797 DO jm = nlay_s, 2, -1 798 DO ji = 1, npti 799 jk = jm - 1 800 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 801 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 802 END DO 803 END DO 772 804 ! 773 805 !-------------------------------------------------------------- … … 923 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 924 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN 925 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji, 1) &926 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &957 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s) & 958 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 927 959 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 928 960 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90
r13642 r14026 95 95 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 96 96 !!--------------------------------------------------------------------- 97 IF( ln_timing ) CALL timing_start('ice _update')97 IF( ln_timing ) CALL timing_start('iceupdate') 98 98 99 99 IF( kt == nit000 .AND. lwp ) THEN … … 156 156 ! ice-ocean mass flux 157 157 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 158 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)158 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 159 159 160 160 ! snw-ocean mass flux … … 162 162 163 163 ! total mass flux at the ocean/ice interface 164 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model165 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux164 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 165 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 166 166 167 167 ! Salt flux at the ocean surface … … 174 174 snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step 175 175 ! ! new mass per unit area 176 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )176 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * ( vt_ip(ji,jj) + vt_il(ji,jj) ) ) 177 177 ! ! time evolution of snow+ice mass 178 178 snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice … … 289 289 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 290 290 IF( ln_ctl ) CALL ice_prt3D ('iceupdate') ! prints 291 IF( ln_timing ) CALL timing_stop ('ice _update')! timing291 IF( ln_timing ) CALL timing_stop ('iceupdate') ! timing 292 292 ! 293 293 END SUBROUTINE ice_update_flx … … 327 327 REAL(wp) :: zflagi ! - - 328 328 !!--------------------------------------------------------------------- 329 IF( ln_timing ) CALL timing_start('ice _update_tau')329 IF( ln_timing ) CALL timing_start('iceupdate') 330 330 331 331 IF( kt == nit000 .AND. lwp ) THEN … … 383 383 CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1., vtau, 'V', -1. ) ! lateral boundary condition 384 384 ! 385 IF( ln_timing ) CALL timing_stop('ice _update_tau')385 IF( ln_timing ) CALL timing_stop('iceupdate') 386 386 ! 387 387 END SUBROUTINE ice_update_tau -
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icevar.F90
r13284 r14026 228 228 z1_zhmax = 1._wp / hi_max(jpl) 229 229 WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 230 h_i (:,:,jpl) = zhmax230 h_i (:,:,jpl) = zhmax 231 231 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 232 232 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) … … 244 244 ELSEWHERE( h_il(:,:,:) >= zhl_max ) ; a_ip_eff(:,:,:) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow 245 245 ELSEWHERE ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * & ! lid is in between. Expose part of the pond 246 & ( h_il(:,:,:) - zhl_min) / ( zhl_max - zhl_min )246 & ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 247 247 END WHERE 248 248 ! … … 545 545 DO ji = 1 , jpi 546 546 ! update exchanges with ocean 547 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice 548 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice 549 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice 547 sfx_res(ji,jj) = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice 548 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice 549 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice 550 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_rdtice 550 551 ! 551 552 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) … … 562 563 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 563 564 v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 565 h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 566 h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 564 567 ! 565 568 END DO … … 656 659 psv_i (ji,jj,jl) = 0._wp 657 660 ENDIF 661 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 662 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 663 pv_il (ji,jj,jl) = 0._wp 664 ENDIF 665 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 666 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 667 pv_ip (ji,jj,jl) = 0._wp 668 ENDIF 658 669 END DO 659 670 END DO … … 665 676 WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp 666 677 WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp 667 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)668 WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok669 678 ! 670 679 END SUBROUTINE ice_var_zapneg
Note: See TracChangeset
for help on using the changeset viewer.