- Timestamp:
- 2015-02-09T14:39:07+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5067 r5070 160 160 zremap_flag(ji,jj) = 0 161 161 ENDIF 162 END DO !ji163 END DO !jj162 END DO 163 END DO 164 164 165 165 !----------------------------------------------------------------------------------------------- … … 201 201 END DO 202 202 203 END DO !jl203 END DO 204 204 205 205 !----------------------------------------------------------------------------------------------- … … 244 244 !----------------------------------------------------------------------------------------------- 245 245 !- 7.1 g(h) for category 1 at start of time step 246 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 247 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 246 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 248 247 & hR(:,:,klbnd), zremap_flag ) 249 248 … … 254 253 255 254 !ji 256 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN255 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 257 256 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 258 257 ! ji, a_i > epsi10 259 IF (zdh0 .lt. 0.0) THEN !remove area from category 1258 IF( zdh0 < 0.0 ) THEN !remove area from category 1 260 259 ! ji, a_i > epsi10; zdh0 < 0 261 zdh0 = MIN( -zdh0,hi_max(klbnd))260 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 262 261 263 262 !Integrate g(1) from 0 to dh0 to estimate area melted 264 zetamax = MIN( zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd)265 IF (zetamax.gt.0.0) THEN263 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 264 IF( zetamax > 0.0 ) THEN 266 265 zx1 = zetamax 267 zx2 = 0.5 * zetamax *zetamax266 zx2 = 0.5 * zetamax * zetamax 268 267 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 269 268 ! Constrain new thickness <= ht_i 270 zdamax = a_i(ii,ij,klbnd) * & 271 (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 272 270 !ice area lost due to melting of thin ice 273 zda0 = MIN( zda0, zdamax)271 zda0 = MIN( zda0, zdamax ) 274 272 275 273 ! Remove area, conserving volume 276 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 277 * 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 ) 278 275 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 279 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) *ht_i(ii,ij,klbnd) ! clem-useless ?280 ENDIF ! zetamax > 0276 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 277 ENDIF 281 278 ! ji, a_i > epsi10 282 279 283 280 ELSE ! if ice accretion 284 281 ! ji, a_i > epsi10; zdh0 > 0 285 zhbnew(ii,ij,klbnd-1) = MIN( zdh0,hi_max(klbnd))282 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 286 283 ! zhbnew was 0, and is shifted to the right to account for thin ice 287 284 ! growth in openwater (F0 = f1) … … 295 292 !- 7.3 g(h) for each thickness category 296 293 DO jl = klbnd, kubnd 297 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&298 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 ) 299 296 END DO 300 297 … … 316 313 ij = nind_j(ji) 317 314 318 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 319 316 320 317 ! left and right integration limits in eta space 321 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)322 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) 323 320 zdonor(ii,ij,jl) = jl 324 321 … … 327 324 ! left and right integration limits in eta space 328 325 zvetamin(ji) = 0.0 329 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) 330 327 zdonor(ii,ij,jl) = jl + 1 331 328 332 329 ENDIF ! zhbnew(jl) > hi_max(jl) 333 330 334 zetamax = MAX( zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin331 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 335 332 zetamin = zvetamin(ji) 336 333 337 334 zx1 = zetamax - zetamin 338 zwk1 = zetamin *zetamin339 zwk2 = zetamax *zetamax340 zx2 = 0.5 * ( zwk2 - zwk1)335 zwk1 = zetamin * zetamin 336 zwk2 = zetamax * zetamax 337 zx2 = 0.5 * ( zwk2 - zwk1 ) 341 338 zwk1 = zwk1 * zetamin 342 339 zwk2 = zwk2 * zetamax 343 zx3 = 1.0 /3.0 * (zwk2 - zwk1)340 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 344 341 nd = zdonor(ii,ij,jl) 345 342 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 346 343 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 347 344 348 END DO ! ji345 END DO 349 346 END DO ! jl klbnd -> kubnd - 1 350 347 … … 365 362 ht_i(ii,ij,1) = rn_himin 366 363 ENDIF 367 END DO !ji364 END DO 368 365 369 366 !!---------------------------------------------------------------------------------------------- … … 401 398 402 399 403 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 404 & g0, g1, hL, hR, zremap_flag ) 400 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 405 401 !!------------------------------------------------------------------ 406 402 !! *** ROUTINE lim_itd_fitline *** … … 442 438 ! Change hL or hR if hice falls outside central third of range 443 439 444 zh13 = 1.0 /3.0 * (2.0*hL(ji,jj) + hR(ji,jj))445 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) ) 446 442 447 443 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 454 450 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 455 451 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 456 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )457 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 ) 458 454 ! 459 455 ELSE ! remap_flag = .false. or a_i < epsi10 … … 516 512 517 513 DO jl = klbnd, kubnd 518 zaTsfn(:,:,jl) = a_i(:,:,jl) *t_su(:,:,jl)514 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 519 515 END DO 520 516 … … 539 535 DO ji = 1, jpi 540 536 541 IF (zdonor(ji,jj,jl) .GT.0) THEN537 IF (zdonor(ji,jj,jl) > 0) THEN 542 538 jl1 = zdonor(ji,jj,jl) 543 539 544 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 545 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 546 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 547 .OR. & 548 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 549 ) 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 550 544 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 551 545 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 559 553 ENDIF 560 554 561 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 562 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 563 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 564 .OR. & 565 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 566 ) 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 567 559 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 568 560 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 577 569 578 570 ! If daice is close to aicen, set daice = aicen. 579 IF (zdaice(ji,jj,jl) .GT.a_i(ji,jj,jl1) - epsi10 ) THEN580 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 581 573 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 582 574 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 586 578 ENDIF 587 579 588 IF (zdvice(ji,jj,jl) .GT.v_i(ji,jj,jl1)-epsi10) THEN589 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 590 582 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 591 583 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) … … 596 588 597 589 ENDIF ! donor > 0 598 END DO ! i599 END DO ! j600 601 END DO !jl590 END DO 591 END DO 592 593 END DO 602 594 603 595 !------------------------------------------------------------------------------- … … 609 601 DO jj = 1, jpj 610 602 DO ji = 1, jpi 611 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 612 604 nbrem = nbrem + 1 613 605 nind_i(nbrem) = ji 614 606 nind_j(nbrem) = jj 615 ENDIF ! tmask607 ENDIF 616 608 END DO 617 609 END DO … … 622 614 623 615 jl1 = zdonor(ii,ij,jl) 624 rswitch 625 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 626 618 IF( jl1 == jl) THEN ; jl2 = jl1+1 627 ELSE 619 ELSE ; jl2 = jl 628 620 ENDIF 629 621 … … 682 674 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 683 675 684 END DO ! ji676 END DO 685 677 686 678 !------------------ … … 689 681 690 682 DO jk = 1, nlay_i 691 !CDIR NODEP692 683 DO ji = 1, nbrem 693 684 ii = nind_i(ji) … … 695 686 696 687 jl1 = zdonor(ii,ij,jl) 697 IF (jl1 .EQ.jl) THEN688 IF (jl1 == jl) THEN 698 689 jl2 = jl+1 699 690 ELSE ! n1 = n+1 … … 704 695 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 705 696 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 706 END DO ! ji707 END DO ! jk697 END DO 698 END DO 708 699 709 700 END DO ! boundaries, 1 to ncat-1 … … 719 710 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 720 711 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 721 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 722 713 ELSE 723 714 ht_i(ji,jj,jl) = 0._wp 724 715 t_su(ji,jj,jl) = rtt 725 716 ENDIF 726 END DO ! ji727 END DO ! jj728 END DO ! jl717 END DO 718 END DO 719 END DO 729 720 ! 730 721 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) … … 836 827 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 837 828 ENDIF 838 END DO ! ji839 END DO ! jj829 END DO 830 END DO 840 831 IF(lk_mpp) CALL mpp_max( zshiftflag ) 841 832 … … 870 861 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 871 862 ENDIF 872 END DO ! ji873 END DO ! jj863 END DO 864 END DO 874 865 875 866 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 890 881 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 891 882 ! ENDIF 892 ! END DO ! ji893 ! END DO ! jj883 ! END DO 884 ! END DO 894 885 ! clem-change end 895 886 896 END DO ! jl887 END DO 897 888 898 889 !------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
