- Timestamp:
- 2017-09-26T15:24:17+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
r8534 r8563 101 101 zhbnew(:,:) = 0._wp 102 102 103 CALL tab_3d_2d( nidx, idxice(1:nidx), h t_i_2d (1:nidx,1:jpl), ht_i )104 CALL tab_3d_2d( nidx, idxice(1:nidx), h t_ib_2d(1:nidx,1:jpl), ht_i_b )103 CALL tab_3d_2d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i ) 104 CALL tab_3d_2d( nidx, idxice(1:nidx), h_ib_2d(1:nidx,1:jpl), h_i_b ) 105 105 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 106 106 CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b ) … … 109 109 ! Compute thickness change in each ice category 110 110 DO ji = 1, nidx 111 zdhice(ji,jl) = h t_i_2d(ji,jl) - ht_ib_2d(ji,jl)111 zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl) 112 112 END DO 113 113 END DO … … 119 119 ! 120 120 ! --- New boundary: Hn* = Hn + Fn*dt --- ! 121 ! Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h t_i_b)121 ! Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h_i_b) 122 122 ! 123 123 IF ( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl+1) & a(jl) /= 0 124 zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h t_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) )125 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h t_ib_2d(ji,jl) )124 zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) ) 125 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) ) 126 126 ELSEIF( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN ! a(jl+1)=0 => Hn* = Hn + fn*dt 127 127 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) … … 136 136 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 137 137 ! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 138 IF( a_i_2d(ji,jl ) > epsi10 .AND. h t_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) idxice(ji) = 0139 IF( a_i_2d(ji,jl+1) > epsi10 .AND. h t_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) idxice(ji) = 0138 IF( a_i_2d(ji,jl ) > epsi10 .AND. h_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) idxice(ji) = 0 139 IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) idxice(ji) = 0 140 140 141 141 ! 2) Hn-1 < Hn* < Hn+1 … … 149 149 DO ji = 1, nidx 150 150 IF( a_i_2d(ji,jpl) > epsi10 ) THEN 151 zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h t_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )151 zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) 152 152 ELSE 153 153 zhbnew(ji,jpl) = hi_max(jpl) … … 158 158 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 159 159 ! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 160 IF( h t_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) idxice(ji) = 0161 IF( h t_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) idxice(ji) = 0160 IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) idxice(ji) = 0 161 IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) idxice(ji) = 0 162 162 END DO 163 163 ! … … 189 189 DO jl = 1, jpl 190 190 ! 191 CALL tab_2d_1d( nidx, idxice(1:nidx), h t_ib_1d(1:nidx), ht_i_b(:,:,jl) )192 CALL tab_2d_1d( nidx, idxice(1:nidx), h t_i_1d (1:nidx), ht_i(:,:,jl) )191 CALL tab_2d_1d( nidx, idxice(1:nidx), h_ib_1d(1:nidx), h_i_b(:,:,jl) ) 192 CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) ) 193 193 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 194 194 CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) … … 197 197 ! 198 198 ! --- g(h) for category 1 --- ! 199 CALL ice_itd_glinear( zhb0(1:nidx) , zhb1(1:nidx) , h t_ib_1d(1:nidx) , a_i_1d(1:nidx) , & ! in199 CALL ice_itd_glinear( zhb0(1:nidx) , zhb1(1:nidx) , h_ib_1d(1:nidx) , a_i_1d(1:nidx) , & ! in 200 200 & g0 (1:nidx,1), g1 (1:nidx,1), hL (1:nidx,1), hR (1:nidx,1) ) ! out 201 201 ! … … 205 205 IF( a_i_1d(ji) > epsi10 ) THEN 206 206 ! 207 zdh0 = h t_i_1d(ji) - ht_ib_1d(ji)207 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 208 208 IF( zdh0 < 0.0 ) THEN !remove area from category 1 209 209 zdh0 = MIN( -zdh0, hi_max(1) ) … … 215 215 zx2 = 0.5 * zetamax * zetamax 216 216 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 217 zdamax = a_i_1d(ji) * (1.0 - h t_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i217 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 218 218 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 219 219 ! of thin ice (zdamax > 0) 220 220 ! Remove area, conserving volume 221 h t_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )221 h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 222 222 a_i_1d(ji) = a_i_1d(ji) - zda0 223 v_i_1d(ji) = a_i_1d(ji) * h t_i_1d(ji) ! clem-useless ?223 v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! clem-useless ? 224 224 ENDIF 225 225 ! … … 233 233 END DO 234 234 ! 235 CALL tab_1d_2d( nidx, idxice(1:nidx), h t_i_1d (1:nidx), ht_i(:,:,jl) )235 CALL tab_1d_2d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) ) 236 236 CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 237 237 CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) … … 240 240 ! 241 241 ! --- g(h) for each thickness category --- ! 242 CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), h t_i_1d(1:nidx) , a_i_1d(1:nidx) , & ! in242 CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), h_i_1d(1:nidx) , a_i_1d(1:nidx) , & ! in 243 243 & g0 (1:nidx,jl ), g1 (1:nidx,jl), hL (1:nidx,jl), hR (1:nidx,jl) ) ! out 244 244 ! … … 284 284 285 285 !---------------------------------------------------------------------------------------------- 286 ! 7) Make sure h t_i >= minimum ice thickness hi_min286 ! 7) Make sure h_i >= minimum ice thickness hi_min 287 287 !---------------------------------------------------------------------------------------------- 288 CALL tab_2d_1d( nidx, idxice(1:nidx), h t_i_1d (1:nidx), ht_i(:,:,1) )288 CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1) ) 289 289 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) ) 290 290 CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) ) 291 291 292 292 DO ji = 1, nidx 293 IF ( a_i_1d(ji) > epsi10 .AND. h t_i_1d(ji) < rn_himin ) THEN294 a_i_1d (ji) = a_i_1d(ji) * h t_i_1d(ji) / rn_himin293 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 294 a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 295 295 ! MV MP 2016 296 296 IF ( nn_pnd_scheme > 0 ) THEN 297 a_ip_1d(ji) = a_ip_1d(ji) * h t_i_1d(ji) / rn_himin297 a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 298 298 ENDIF 299 299 ! END MV MP 2016 300 h t_i_1d(ji) = rn_himin300 h_i_1d(ji) = rn_himin 301 301 ENDIF 302 302 END DO 303 303 ! 304 CALL tab_1d_2d( nidx, idxice(1:nidx), h t_i_1d (1:nidx), ht_i(:,:,1) )304 CALL tab_1d_2d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1) ) 305 305 CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) ) 306 306 CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) ) … … 392 392 !!------------------------------------------------------------------ 393 393 394 CALL tab_3d_2d( nidx, idxice(1:nidx), h t_i_2d (1:nidx,1:jpl), ht_i )394 CALL tab_3d_2d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i ) 395 395 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 396 396 CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i ) … … 518 518 !------------------------------------------------------------------------------- 519 519 WHERE( a_i_2d(1:nidx,:) >= epsi20 ) 520 h t_i_2d(1:nidx,:) = v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:)520 h_i_2d(1:nidx,:) = v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 521 521 t_su_2d(1:nidx,:) = zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 522 522 ELSEWHERE 523 h t_i_2d(1:nidx,:) = 0._wp523 h_i_2d(1:nidx,:) = 0._wp 524 524 t_su_2d(1:nidx,:) = rt0 525 525 END WHERE 526 526 ! 527 CALL tab_2d_3d( nidx, idxice(1:nidx), h t_i_2d (1:nidx,1:jpl), ht_i )527 CALL tab_2d_3d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i ) 528 528 CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 529 529 CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i ) … … 574 574 END DO 575 575 ! 576 !!clem CALL tab_2d_1d( nidx, idxice(1:nidx), h t_i_1d (1:nidx), ht_i(:,:,jl) )576 !!clem CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl) ) 577 577 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 578 578 CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) … … 582 582 ! how much of a_i you send in cat sup is somewhat arbitrary 583 583 !!clem: these do not work properly after a restart (I do not know why) 584 !! zdaice(ji,jl) = a_i_1d(ji) * ( h t_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji)584 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 585 585 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 586 586 !!clem: these do not work properly after a restart (I do not know why)
Note: See TracChangeset
for help on using the changeset viewer.
