Changeset 868 for trunk/NEMO/LIM_SRC_3/limrhg.F90
- Timestamp:
- 2008-03-14T19:53:00+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limrhg.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limrhg.F90
r866 r868 31 31 PUBLIC lim_rhg ! routine called by lim_dyn 32 32 33 !! * Substitutions 34 # include "vectopt_loop_substitute.h90" 35 33 36 !! * Module variables 34 37 REAL(wp) :: & ! constant values … … 111 114 112 115 INTEGER :: & 113 iter, jter!: temporary integers116 jter !: temporary integers 114 117 115 118 CHARACTER (len=50) :: charout … … 182 185 ! 183 186 ! Put every vector to 0 184 zpresh(:,:) = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:) = 0.0 185 zmass1(:,:) = 0.0 ; zcorl1(:,:) = 0.0 ; zcorl2(:,:) = 0.0 186 za1ct(:,:) = 0.0 ; za2ct(:,:) = 0.0 187 zc1(:,:) = 0.0 ; zusw(:,:) = 0.0 188 u_oce1(:,:) = 0.0 ; v_oce1(:,:) = 0.0 ; u_oce2(:,:) = 0.0 189 v_oce2(:,:) = 0.0 ; u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0 190 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 187 zpresh(:,:) = 0.0 ; zc1(:,:) = 0.0 188 zpreshc(:,:) = 0.0 189 u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0 191 190 zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0 192 deltat(:,:) = 0.0 ; deltac(:,:) = 0.0193 zpresh(:,:) = 0.0194 191 195 192 ! Ice strength on T-points … … 197 194 198 195 ! Ice mass and temp variables 199 DO jj = k_j1 , k_jpj-1 196 !CDIR NOVERRCHK 197 DO jj = k_j1 , k_jpj 198 !CDIR NOVERRCHK 200 199 DO ji = 1 , jpi 201 200 zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) … … 207 206 END DO 208 207 209 CALL lbc_lnk( zc1(:,:), 'T', 1. )210 CALL lbc_lnk( zpresh(:,:), 'T', 1. )211 212 CALL lbc_lnk( tmi(:,:), 'T', 1. )213 214 208 ! Ice strength on grid cell corners (zpreshc) 215 209 ! needed for calculation of shear stress 210 !CDIR NOVERRCHK 216 211 DO jj = k_j1+1, k_jpj-1 217 DO ji = 2, jpim1 212 !CDIR NOVERRCHK 213 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 218 214 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 219 215 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & … … 249 245 250 246 DO jj = k_j1+1, k_jpj-1 251 DO ji = 2,jpim1247 DO ji = fs_2, fs_jpim1 252 248 253 249 zt11 = tms(ji,jj)*e1t(ji,jj) … … 276 272 277 273 ! Ocean has no slip boundary condition 278 ! GG bug279 ! v_oce1(ji,jj) = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj) &280 ! & +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) &281 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)282 274 v_oce1(ji,jj) = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji,jj) & 283 275 & +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji+1,jj)) & 284 276 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 285 277 286 ! GG bug287 ! u_oce2(ji,jj) = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1) &288 ! & +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) &289 ! & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)290 278 u_oce2(ji,jj) = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj) & 291 279 & +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj+1)) & … … 330 318 zs12(:,:) = stress12_i(:,:) 331 319 332 v_ice1(:,:) = 0.0333 u_ice2(:,:) = 0.0334 335 zdd(:,:) = 0.0336 zdt(:,:) = 0.0337 zds(:,:) = 0.0338 deltat(:,:) = 0.0339 320 !----------------------! 340 DO iter = 1 , nevp ! loop over iter !321 DO jter = 1 , nevp ! loop over jter ! 341 322 !----------------------! 342 323 DO jj = k_j1, k_jpj-1 … … 346 327 347 328 DO jj = k_j1+1, k_jpj-1 348 DO ji = 2,jpim1329 DO ji = fs_2, fs_jpim1 349 330 350 331 ! … … 407 388 END DO 408 389 409 CALL lbc_lnk( zdd(:,:), 'T', 1. ) 410 CALL lbc_lnk( zdt(:,:), 'T', 1. ) 411 CALL lbc_lnk( zds(:,:), 'F', 1. ) 412 390 !CDIR NOVERRCHK 413 391 DO jj = k_j1+1, k_jpj-1 414 DO ji = 2, jpim1 392 !CDIR NOVERRCHK 393 DO ji = fs_2, fs_jpim1 415 394 416 395 !- Calculate Delta at centre of grid cells … … 446 425 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 447 426 427 !CDIR NOVERRCHK 448 428 DO jj = k_j1+1, k_jpj-1 449 DO ji = 2, jpim1 429 !CDIR NOVERRCHK 430 DO ji = fs_2, fs_jpim1 450 431 !- Calculate Delta on corners 451 432 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & … … 482 463 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 483 464 DO jj = k_j1+1, k_jpj-1 484 DO ji = 2,jpim1465 DO ji = fs_2, fs_jpim1 485 466 !- contribution of zs1, zs2 and zs12 to zf1 486 467 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & … … 501 482 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 502 483 ! 503 IF (MOD(iter,2).eq.0) THEN 504 484 IF (MOD(jter,2).eq.0) THEN 485 486 !CDIR NOVERRCHK 505 487 DO jj = k_j1+1, k_jpj-1 506 DO ji = 2, jpim1 488 !CDIR NOVERRCHK 489 DO ji = fs_2, fs_jpim1 507 490 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 508 491 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg … … 510 493 511 494 ! SB modif because ocean has no slip boundary condition 512 ! GG bug513 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) &514 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &515 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)516 495 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 517 496 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & … … 530 509 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 531 510 511 !CDIR NOVERRCHK 532 512 DO jj = k_j1+1, k_jpj-1 533 DO ji = 2, jpim1 534 535 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 536 zsang = sign(1.0,fcor(ji,jj))*sangvg 513 !CDIR NOVERRCHK 514 DO ji = fs_2, fs_jpim1 515 516 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 517 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 537 518 z0 = zmass2(ji,jj)/dtevp 538 519 ! SB modif because ocean has no slip boundary condition 539 ! GG bug540 ! zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) &541 ! & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &542 ! & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)543 520 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 544 521 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & … … 558 535 559 536 ELSE 537 !CDIR NOVERRCHK 560 538 DO jj = k_j1+1, k_jpj-1 561 DO ji = 2, jpim1 562 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 563 zsang = sign(1.0,fcor(ji,jj))*sangvg 539 !CDIR NOVERRCHK 540 DO ji = fs_2, fs_jpim1 541 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 542 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 564 543 z0 = zmass2(ji,jj)/dtevp 565 544 ! SB modif because ocean has no slip boundary condition 566 ! GG Bug567 ! zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) &568 ! & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &569 ! & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)570 545 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 571 546 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & … … 585 560 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 586 561 562 !CDIR NOVERRCHK 587 563 DO jj = k_j1+1, k_jpj-1 588 DO ji = 2, jpim1 564 !CDIR NOVERRCHK 565 DO ji = fs_2, fs_jpim1 589 566 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 590 567 zsang = SIGN(1.0,fcor(ji,jj))*sangvg … … 613 590 ENDIF 614 591 615 !--- Convergence test. 616 DO jj = k_j1+1 , k_jpj-1 617 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 618 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 619 END DO 620 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 592 IF(ln_ctl) THEN 593 !--- Convergence test. 594 DO jj = k_j1+1 , k_jpj-1 595 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 596 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 597 END DO 598 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 599 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 600 ENDIF 621 601 622 602 ! ! ==================== ! 623 END DO ! end loop over iter !603 END DO ! end loop over jter ! 624 604 ! ! ==================== ! 625 605 … … 632 612 ! ocean velocity, 633 613 ! This prevents high velocity when ice is thin 614 !CDIR NOVERRCHK 634 615 DO jj = k_j1+1, k_jpj-1 635 DO ji = 2, jpim1 616 !CDIR NOVERRCHK 617 DO ji = fs_2, fs_jpim1 636 618 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 637 619 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) … … 643 625 END DO 644 626 627 DO jj = k_j1+1, k_jpj-1 628 DO ji = fs_2, fs_jpim1 629 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 630 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 631 IF ( zdummy .LE. 5.0e-2 ) THEN 632 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 633 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 634 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 635 636 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 637 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 638 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 639 ENDIF ! zdummy 640 END DO 641 END DO 642 643 CALL lbc_lnk( u_ice2(:,:), 'U', -1. ) 644 CALL lbc_lnk( v_ice1(:,:), 'V', -1. ) 645 645 646 ! Recompute delta, shear and div, inputs for mechanical redistribution 647 !CDIR NOVERRCHK 646 648 DO jj = k_j1+1, k_jpj-1 647 DO ji = 2, jpim1 649 !CDIR NOVERRCHK 650 DO ji = fs_2, fs_jpim1 648 651 !- zdd(:,:), zdt(:,:): divergence and tension at centre 649 652 !- zds(:,:): shear on northeast corner of grid cells … … 680 683 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 681 684 682 !- Calculate Delta at centre of grid cells683 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) &684 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &685 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)686 687 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) &688 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &689 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)690 691 685 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & 692 686 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & … … 710 704 ! 711 705 ! * Invariants of the stress tensor are required for limitd_me 712 ! * Store the stress tensor for the next time step713 706 ! accelerates convergence and improves stability 714 707 DO jj = k_j1+1, k_jpj-1 715 DO ji = 2, jpim1 716 divu_i(ji,jj) = zdd(ji,jj) 717 delta_i(ji,jj) = deltat (ji,jj) 718 shear_i(ji,jj) = zds (ji,jj) 719 stress1_i(ji,jj) = zs1(ji,jj) 720 stress2_i(ji,jj) = zs2(ji,jj) 721 stress12_i(ji,jj) = zs12(ji,jj) 708 DO ji = fs_2, fs_jpim1 709 divu_i (ji,jj) = zdd (ji,jj) 710 delta_i(ji,jj) = deltat(ji,jj) 711 shear_i(ji,jj) = zds (ji,jj) 722 712 END DO 723 713 END DO 724 714 725 715 ! Lateral boundary condition 726 CALL lbc_lnk( divu_i (:,:), 'T', 1. )716 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 727 717 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 728 718 CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 729 CALL lbc_lnk( stress1_i(:,:), 'T', 1. ) 730 CALL lbc_lnk( stress2_i(:,:), 'T', 1. ) 731 CALL lbc_lnk( stress12_i(:,:), 'F', 1. ) 719 720 ! * Store the stress tensor for the next time step 721 stress1_i (:,:) = zs1 (:,:) 722 stress2_i (:,:) = zs2 (:,:) 723 stress12_i(:,:) = zs12(:,:) 732 724 733 725 ! … … 738 730 ! print the residual for convergence 739 731 IF(ln_ctl) THEN 740 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, iter732 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 741 733 CALL prt_ctl_info(charout) 742 734 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
Note: See TracChangeset
for help on using the changeset viewer.
