Changeset 9880 for NEMO/trunk/src/ICE
- Timestamp:
- 2018-07-05T17:43:44+02:00 (6 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 3 edited
-
icedyn_adv.F90 (modified) (1 diff)
-
icedyn_rdgrft.F90 (modified) (1 diff)
-
iceitd.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv.F90
r9604 r9880 95 95 END SELECT 96 96 97 !---------------------------- 98 ! Debug the advection schemes 99 !---------------------------- 100 ! clem: The 2 advection schemes above are not strictly positive. 101 ! In Prather, advected fields are bounded by 0 in the routine with a MAX(0,field) ==> likely conservation issues 102 ! In UMx , advected fields are not bounded and negative values can appear. 103 ! These values are usually very small but in some occasions they can also be non-negligible 104 ! Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 105 ! ==> 1) remove negative ice areas and volumes (conservation is ensure) 106 CALL ice_var_zapsmall 107 ! ==> 2) remove remaining negative advected fields (conservation is not preserved) 108 WHERE( v_s (:,:,:) < 0._wp ) v_s (:,:,:) = 0._wp 109 WHERE( sv_i(:,:,:) < 0._wp ) sv_i(:,:,:) = 0._wp 110 WHERE( e_i (:,:,:,:) < 0._wp ) e_i (:,:,:,:) = 0._wp 111 WHERE( e_s (:,:,:,:) < 0._wp ) e_s (:,:,:,:) = 0._wp 112 97 113 !------------ 98 114 ! diagnostics -
NEMO/trunk/src/ICE/icedyn_rdgrft.F90
r9604 r9880 143 143 INTEGER, PARAMETER :: jp_itermax = 20 144 144 !!------------------------------------------------------------------- 145 ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 146 ! likely due to truncation error ( i.e. 1. - 1. /= 0 ) 147 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 148 145 149 ! controls 146 150 IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing -
NEMO/trunk/src/ICE/iceitd.F90
r9604 r9880 95 95 DO ji = 1, jpi 96 96 IF ( at_i(ji,jj) > epsi10 ) THEN 97 npti 97 npti = npti + 1 98 98 nptidx( npti ) = (jj - 1) * jpi + ji 99 99 ENDIF … … 111 111 CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i ) 112 112 CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b ) 113 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i)114 CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b)113 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 114 CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d(1:npti,1:jpl), a_i_b ) 115 115 ! 116 116 DO jl = 1, jpl 117 117 ! Compute thickness change in each ice category 118 118 DO ji = 1, npti 119 IF( a_i_2d(ji,jl) > epsi10 ) zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl) ! clem: if statement is necessary for repro119 IF( a_i_2d(ji,jl) > epsi10 ) zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl) 120 120 END DO 121 121 END DO … … 130 130 ! 131 131 IF ( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl+1) & a(jl) /= 0 132 zslope 132 zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) ) 133 133 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) ) 134 134 ELSEIF( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN ! a(jl+1)=0 => Hn* = Hn + fn*dt … … 198 198 ! 199 199 CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) ) 200 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl))201 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl))202 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl))200 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 201 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 202 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,jl) ) 203 203 ! 204 204 IF( jl == 1 ) THEN … … 206 206 ! --- g(h) for category 1 --- ! 207 207 CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in 208 & g0 (1:npti,1), g1 (1:npti,1), hL 208 & g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out 209 209 ! 210 210 ! Area lost due to melting of thin ice … … 228 228 ! Remove area, conserving volume 229 229 h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 230 a_i_1d(ji) 231 v_i_1d(ji) 230 a_i_1d(ji) = a_i_1d(ji) - zda0 231 v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! useless ? 232 232 ENDIF 233 233 ! … … 241 241 END DO 242 242 ! 243 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl))244 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl))245 CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl))243 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 244 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 245 CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 246 246 ! 247 247 ENDIF ! jl=1 … … 249 249 ! --- g(h) for each thickness category --- ! 250 250 CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in 251 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR 251 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out 252 252 ! 253 253 END DO … … 300 300 DO ji = 1, npti 301 301 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 302 a_i_1d 302 a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 303 303 IF( ln_pnd_H12 ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 304 304 h_i_1d(ji) = rn_himin … … 436 436 ENDIF 437 437 ! 438 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 439 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 440 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 441 ! 438 442 a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas 439 443 a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl) … … 445 449 v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 446 450 v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 447 ! 448 ! ! Ice age 449 ztrans = oa_i_2d(ji,jl1) * zworka(ji) 451 ! 452 ztrans = oa_i_2d(ji,jl1) * zworka(ji) ! Ice age 450 453 oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans 451 454 oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans … … 455 458 sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans 456 459 ! 457 ! ! Surface temperature 458 ztrans = zaTsfn(ji,jl1) * zworka(ji) 460 ztrans = zaTsfn(ji,jl1) * zworka(ji) ! Surface temperature 459 461 zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans 460 462 zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans 461 463 ! 462 464 IF ( ln_pnd_H12 ) THEN 463 ! ! Pond fraction 464 ztrans = a_ip_2d(ji,jl1) * zworka(ji) 465 ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction 465 466 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 466 467 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 467 ! ! Pond volume (also proportional to da/a)468 ztrans = v_ip_2d(ji,jl1) * zworka(ji) 468 ! 469 ztrans = v_ip_2d(ji,jl1) * zworka(ji) ! Pond volume (also proportional to da/a) 469 470 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 470 471 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans … … 519 520 !------------------------------------------------------------------------------- 520 521 WHERE( a_i_2d(1:npti,:) >= epsi20 ) 521 h_i_2d (1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:)522 h_i_2d (1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 522 523 t_su_2d(1:npti,:) = zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 523 524 ELSEWHERE 524 h_i_2d (1:npti,:) = 0._wp525 h_i_2d (1:npti,:) = 0._wp 525 526 t_su_2d(1:npti,:) = rt0 526 527 END WHERE … … 568 569 DO jj = 1, jpj 569 570 DO ji = 1, jpi 570 IF( a_i(ji,jj,jl) > epsi10.AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN571 IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 571 572 npti = npti + 1 572 573 nptidx( npti ) = (jj - 1) * jpi + ji … … 575 576 END DO 576 577 ! 577 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl))578 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl))579 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl))578 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 579 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 580 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 580 581 ! 581 582 DO ji = 1, npti 582 583 jdonor(ji,jl) = jl 583 584 ! how much of a_i you send in cat sup is somewhat arbitrary 584 !!clem: these do not work properly after a restart (I do not know why) 585 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 585 586 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 586 587 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 587 !!clem: these do not work properly after a restart (I do not know why) 588 !! 589 !! 588 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 589 !! zdaice(ji,jl) = a_i_1d(ji) 590 !! zdvice(ji,jl) = v_i_1d(ji) 590 591 !!clem: these are from UCL and work ok 591 592 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp … … 609 610 DO jj = 1, jpj 610 611 DO ji = 1, jpi 611 IF( a_i(ji,jj,jl+1) > epsi10.AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN612 IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 612 613 npti = npti + 1 613 614 nptidx( npti ) = (jj - 1) * jpi + ji … … 616 617 END DO 617 618 ! 618 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i(:,:,jl+1)) ! jl+1 is ok619 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i(:,:,jl+1)) ! jl+1 is ok619 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 620 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 620 621 DO ji = 1, npti 621 622 jdonor(ji,jl) = jl + 1
Note: See TracChangeset
for help on using the changeset viewer.
