Changeset 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2015-03-04T17:06:03+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4990 r5123 13 13 !! 'key_lim3' : LIM3 sea-ice model 14 14 !!---------------------------------------------------------------------- 15 !! lim_itd_th : thermodynamics of ice thickness distribution16 15 !! lim_itd_th_rem : 17 16 !! lim_itd_th_reb : … … 25 24 USE thd_ice ! LIM-3 thermodynamic variables 26 25 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters28 USE limthd_lac ! LIM-3 lateral accretion29 26 USE limvar ! LIM-3 variables 30 27 USE limcons ! LIM-3 conservation … … 34 31 USE wrk_nemo ! work arrays 35 32 USE lib_fortran ! to use key_nosignedzero 36 USE timing ! Timing 37 USE limcons ! conservation tests 33 USE limcons ! conservation tests 38 34 39 35 IMPLICIT NONE 40 36 PRIVATE 41 37 42 PUBLIC lim_itd_th ! called by ice_stp43 38 PUBLIC lim_itd_th_rem 44 39 PUBLIC lim_itd_th_reb … … 52 47 !!---------------------------------------------------------------------- 53 48 CONTAINS 54 55 SUBROUTINE lim_itd_th( kt )56 !!------------------------------------------------------------------57 !! *** ROUTINE lim_itd_th ***58 !!59 !! ** Purpose : computes the thermodynamics of ice thickness distribution60 !!61 !! ** Method :62 !!------------------------------------------------------------------63 INTEGER, INTENT(in) :: kt ! time step index64 !65 INTEGER :: ji, jj, jk, jl ! dummy loop index66 !67 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b68 !!------------------------------------------------------------------69 IF( nn_timing == 1 ) CALL timing_start('limitd_th')70 71 ! conservation test72 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)73 74 IF( kt == nit000 .AND. lwp ) THEN75 WRITE(numout,*)76 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'77 WRITE(numout,*) '~~~~~~~~~~~'78 ENDIF79 80 !------------------------------------------------------------------------------|81 ! 1) Transport of ice between thickness categories. |82 !------------------------------------------------------------------------------|83 ! Given thermodynamic growth rates, transport ice between84 ! thickness categories.85 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )86 !87 CALL lim_var_glo2eqv ! only for info88 CALL lim_var_agg(1)89 90 !------------------------------------------------------------------------------|91 ! 3) Add frazil ice growing in leads.92 !------------------------------------------------------------------------------|93 CALL lim_thd_lac94 CALL lim_var_glo2eqv ! only for info95 96 IF(ln_ctl) THEN ! Control print97 CALL prt_ctl_info(' ')98 CALL prt_ctl_info(' - Cell values : ')99 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')100 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th : cell area :')101 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :')102 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :')103 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :')104 DO jl = 1, jpl105 CALL prt_ctl_info(' ')106 CALL prt_ctl_info(' - Category : ', ivar1=jl)107 CALL prt_ctl_info(' ~~~~~~~~~~')108 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ')109 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ')110 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ')111 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ')112 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ')113 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ')114 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ')115 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ')116 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ')117 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ')118 DO jk = 1, nlay_i119 CALL prt_ctl_info(' ')120 CALL prt_ctl_info(' - Layer : ', ivar1=jk)121 CALL prt_ctl_info(' ~~~~~~~')122 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ')123 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ')124 END DO125 END DO126 ENDIF127 !128 ! conservation test129 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)130 !131 IF( nn_timing == 1 ) CALL timing_stop('limitd_th')132 END SUBROUTINE lim_itd_th133 !134 49 135 50 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) … … 153 68 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 154 69 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 155 REAL(wp) :: zx3 , zareamin ! - -70 REAL(wp) :: zx3 156 71 CHARACTER (len = 15) :: fieldid 157 72 … … 188 103 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 189 104 190 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model191 192 105 !!---------------------------------------------------------------------------------------------- 193 106 !! 0) Conservation checkand changes in each ice category … … 216 129 DO jj = 1, jpj 217 130 DO ji = 1, jpi 218 rswitch 131 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 219 132 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 220 rswitch 133 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 221 134 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 222 135 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) … … 239 152 DO jj = 1, jpj 240 153 DO ji = 1, jpi 241 IF ( at_i(ji,jj) .gt. zareamin) THEN154 IF ( at_i(ji,jj) > epsi10 ) THEN 242 155 nbrem = nbrem + 1 243 156 nind_i(nbrem) = ji … … 247 160 zremap_flag(ji,jj) = 0 248 161 ENDIF 249 END DO !ji250 END DO !jj162 END DO 163 END DO 251 164 252 165 !----------------------------------------------------------------------------------------------- … … 254 167 !----------------------------------------------------------------------------------------------- 255 168 !- 4.1 Compute category boundaries 256 ! Tricky trick see limitd_me.F90257 ! will be soon removed, CT258 ! hi_max(kubnd) = 99.259 169 zhbnew(:,:,:) = 0._wp 260 170 … … 291 201 END DO 292 202 293 END DO !jl203 END DO 294 204 295 205 !----------------------------------------------------------------------------------------------- … … 334 244 !----------------------------------------------------------------------------------------------- 335 245 !- 7.1 g(h) for category 1 at start of time step 336 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 337 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 246 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 338 247 & hR(:,:,klbnd), zremap_flag ) 339 248 … … 344 253 345 254 !ji 346 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN255 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 347 256 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 348 257 ! ji, a_i > epsi10 349 IF (zdh0 .lt. 0.0) THEN !remove area from category 1258 IF( zdh0 < 0.0 ) THEN !remove area from category 1 350 259 ! ji, a_i > epsi10; zdh0 < 0 351 zdh0 = MIN( -zdh0,hi_max(klbnd))260 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 352 261 353 262 !Integrate g(1) from 0 to dh0 to estimate area melted 354 zetamax = MIN( zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd)355 IF (zetamax.gt.0.0) THEN263 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 264 IF( zetamax > 0.0 ) THEN 356 265 zx1 = zetamax 357 zx2 = 0.5 * zetamax *zetamax266 zx2 = 0.5 * zetamax * zetamax 358 267 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 359 268 ! Constrain new thickness <= ht_i 360 zdamax = a_i(ii,ij,klbnd) * & 361 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 269 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 362 270 !ice area lost due to melting of thin ice 363 zda0 = MIN( zda0, zdamax)271 zda0 = MIN( zda0, zdamax ) 364 272 365 273 ! Remove area, conserving volume 366 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 367 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 274 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 368 275 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 369 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) *ht_i(ii,ij,klbnd) ! clem-useless ?370 ENDIF ! zetamax > 0276 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 277 ENDIF 371 278 ! ji, a_i > epsi10 372 279 373 280 ELSE ! if ice accretion 374 281 ! ji, a_i > epsi10; zdh0 > 0 375 zhbnew(ii,ij,klbnd-1) = MIN( zdh0,hi_max(klbnd))282 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 376 283 ! zhbnew was 0, and is shifted to the right to account for thin ice 377 284 ! growth in openwater (F0 = f1) … … 385 292 !- 7.3 g(h) for each thickness category 386 293 DO jl = klbnd, kubnd 387 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&388 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag)294 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 295 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 389 296 END DO 390 297 … … 406 313 ij = nind_j(ji) 407 314 408 IF (zhbnew(ii,ij,jl) .gt.hi_max(jl)) THEN ! transfer from jl to jl+1315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 409 316 410 317 ! left and right integration limits in eta space 411 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)412 zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl)318 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 319 zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 413 320 zdonor(ii,ij,jl) = jl 414 321 … … 417 324 ! left and right integration limits in eta space 418 325 zvetamin(ji) = 0.0 419 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1)326 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 420 327 zdonor(ii,ij,jl) = jl + 1 421 328 422 329 ENDIF ! zhbnew(jl) > hi_max(jl) 423 330 424 zetamax = MAX( zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin331 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 425 332 zetamin = zvetamin(ji) 426 333 427 334 zx1 = zetamax - zetamin 428 zwk1 = zetamin *zetamin429 zwk2 = zetamax *zetamax430 zx2 = 0.5 * ( zwk2 - zwk1)335 zwk1 = zetamin * zetamin 336 zwk2 = zetamax * zetamax 337 zx2 = 0.5 * ( zwk2 - zwk1 ) 431 338 zwk1 = zwk1 * zetamin 432 339 zwk2 = zwk2 * zetamax 433 zx3 = 1.0 /3.0 * (zwk2 - zwk1)340 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 434 341 nd = zdonor(ii,ij,jl) 435 342 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 436 343 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 437 344 438 END DO ! ji345 END DO 439 346 END DO ! jl klbnd -> kubnd - 1 440 347 … … 451 358 ii = nind_i(ji) 452 359 ij = nind_j(ji) 453 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim) THEN454 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim455 ht_i(ii,ij,1) = hiclim360 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 361 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 362 ht_i(ii,ij,1) = rn_himin 456 363 ENDIF 457 END DO !ji364 END DO 458 365 459 366 !!---------------------------------------------------------------------------------------------- … … 491 398 492 399 493 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 494 & g0, g1, hL, hR, zremap_flag ) 400 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 495 401 !!------------------------------------------------------------------ 496 402 !! *** ROUTINE lim_itd_fitline *** … … 532 438 ! Change hL or hR if hice falls outside central third of range 533 439 534 zh13 = 1.0 /3.0 * (2.0*hL(ji,jj) + hR(ji,jj))535 zh23 = 1.0 /3.0 * (hL(ji,jj) + 2.0*hR(ji,jj))440 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 441 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 536 442 537 443 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 544 450 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 545 451 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 546 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )547 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5)452 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 453 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 548 454 ! 549 455 ELSE ! remap_flag = .false. or a_i < epsi10 … … 606 512 607 513 DO jl = klbnd, kubnd 608 zaTsfn(:,:,jl) = a_i(:,:,jl) *t_su(:,:,jl)514 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 609 515 END DO 610 516 … … 629 535 DO ji = 1, jpi 630 536 631 IF (zdonor(ji,jj,jl) .GT.0) THEN537 IF (zdonor(ji,jj,jl) > 0) THEN 632 538 jl1 = zdonor(ji,jj,jl) 633 539 634 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 635 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 636 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 637 .OR. & 638 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 639 ) THEN 540 IF (zdaice(ji,jj,jl) < 0.0) THEN 541 IF (zdaice(ji,jj,jl) > -epsi10) THEN 542 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 543 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 640 544 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 641 545 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 649 553 ENDIF 650 554 651 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 652 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 653 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 654 .OR. & 655 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 656 ) THEN 555 IF (zdvice(ji,jj,jl) < 0.0) THEN 556 IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 557 IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. & 558 ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 657 559 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 658 560 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 667 569 668 570 ! If daice is close to aicen, set daice = aicen. 669 IF (zdaice(ji,jj,jl) .GT.a_i(ji,jj,jl1) - epsi10 ) THEN670 IF (zdaice(ji,jj,jl) .LT.a_i(ji,jj,jl1)+epsi10) THEN571 IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 572 IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 671 573 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 672 574 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 676 578 ENDIF 677 579 678 IF (zdvice(ji,jj,jl) .GT.v_i(ji,jj,jl1)-epsi10) THEN679 IF (zdvice(ji,jj,jl) .LT.v_i(ji,jj,jl1)+epsi10) THEN580 IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 581 IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 680 582 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 681 583 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 686 588 687 589 ENDIF ! donor > 0 688 END DO ! i689 END DO ! j690 691 END DO !jl590 END DO 591 END DO 592 593 END DO 692 594 693 595 !------------------------------------------------------------------------------- … … 699 601 DO jj = 1, jpj 700 602 DO ji = 1, jpi 701 IF (zdaice(ji,jj,jl) .GT.0.0 ) THEN ! daice(n) can be < puny603 IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 702 604 nbrem = nbrem + 1 703 605 nind_i(nbrem) = ji 704 606 nind_j(nbrem) = jj 705 ENDIF ! tmask607 ENDIF 706 608 END DO 707 609 END DO … … 712 614 713 615 jl1 = zdonor(ii,ij,jl) 714 rswitch 715 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch616 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 617 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 716 618 IF( jl1 == jl) THEN ; jl2 = jl1+1 717 ELSE 619 ELSE ; jl2 = jl 718 620 ENDIF 719 621 … … 772 674 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 773 675 774 END DO ! ji676 END DO 775 677 776 678 !------------------ … … 779 681 780 682 DO jk = 1, nlay_i 781 !CDIR NODEP782 683 DO ji = 1, nbrem 783 684 ii = nind_i(ji) … … 785 686 786 687 jl1 = zdonor(ii,ij,jl) 787 IF (jl1 .EQ.jl) THEN688 IF (jl1 == jl) THEN 788 689 jl2 = jl+1 789 690 ELSE ! n1 = n+1 … … 794 695 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 795 696 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 796 END DO ! ji797 END DO ! jk697 END DO 698 END DO 798 699 799 700 END DO ! boundaries, 1 to ncat-1 … … 809 710 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 810 711 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 811 rswitch = 1.0 - MAX( 0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes712 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 812 713 ELSE 813 714 ht_i(ji,jj,jl) = 0._wp 814 t_su(ji,jj,jl) = rt t715 t_su(ji,jj,jl) = rt0 815 716 ENDIF 816 END DO ! ji817 END DO ! jj818 END DO ! jl717 END DO 718 END DO 719 END DO 819 720 ! 820 721 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) … … 926 827 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 927 828 ENDIF 928 END DO ! ji929 END DO ! jj829 END DO 830 END DO 930 831 IF(lk_mpp) CALL mpp_max( zshiftflag ) 931 832 … … 951 852 zshiftflag = 0 952 853 953 !clem-change954 854 DO jj = 1, jpj 955 855 DO ji = 1, jpi … … 961 861 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 962 862 ENDIF 963 END DO ! ji964 END DO ! jj863 END DO 864 END DO 965 865 966 866 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 973 873 zdvice(:,:,jl) = 0._wp 974 874 ENDIF 975 !clem-change976 875 977 876 ! ! clem-change begin: why not doing that? … … 982 881 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 983 882 ! ENDIF 984 ! END DO ! ji985 ! END DO ! jj883 ! END DO 884 ! END DO 986 885 ! clem-change end 987 886 988 END DO ! jl887 END DO 989 888 990 889 !------------------------------------------------------------------------------ … … 1013 912 !!---------------------------------------------------------------------- 1014 913 CONTAINS 1015 SUBROUTINE lim_itd_th ! Empty routines1016 END SUBROUTINE lim_itd_th1017 SUBROUTINE lim_itd_th_ini1018 END SUBROUTINE lim_itd_th_ini1019 914 SUBROUTINE lim_itd_th_rem 1020 915 END SUBROUTINE lim_itd_th_rem
Note: See TracChangeset
for help on using the changeset viewer.
