- Timestamp:
- 2015-01-20T15:26:13+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4346 r5038 50 50 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2) 51 51 52 REAL(wp) :: epsi10 = 1.e-10_wp !53 REAL(wp) :: rzero = 0._wp ! constant values54 REAL(wp) :: rone = 1._wp ! constant values55 56 52 !! * Substitutions 57 53 # include "vectopt_loop_substitute.h90" … … 121 117 CHARACTER (len=50) :: charout 122 118 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 123 REAL(wp) :: za, zstms, zsang, zmask ! local scalars 119 REAL(wp) :: za, zstms, zmask ! local scalars 120 REAL(wp) :: zc1, zc2, zc3 ! ice mass 124 121 125 122 REAL(wp) :: dtevp ! time step for subcycling … … 127 124 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 128 125 REAL(wp) :: zu_ice2, zv_ice1 ! 129 REAL(wp) :: zddc, zdtc, zzdst ! delta on corners and on centre 126 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 127 REAL(wp) :: zdst ! shear at the center of the grid point 130 128 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 131 129 REAL(wp) :: sigma1, sigma2 ! internal ice stress 132 130 133 131 REAL(wp) :: zresm ! Maximal error on ice velocity 134 REAL(wp) :: zindb ! ice (1) or not (0)135 132 REAL(wp) :: zdummy ! dummy argument 136 133 REAL(wp) :: zintb, zintn ! dummy argument … … 142 139 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 143 140 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 144 REAL(wp), POINTER, DIMENSION(:,:) :: zc1 ! ice mass145 REAL(wp), POINTER, DIMENSION(:,:) :: zusw ! temporary weight for ice strength computation146 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1 ! ocean u/v component on U points 147 142 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2, v_oce2 ! ocean u/v component on V points … … 149 144 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 150 145 151 REAL(wp), POINTER, DIMENSION(:,:) :: zd d , zdt ! Divergence andtension at centre of grid cells146 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells 152 147 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells 153 REAL(wp), POINTER, DIMENSION(:,:) :: zdst ! Shear on centre of grid cells154 REAL(wp), POINTER, DIMENSION(:,:) :: deltat, deltac ! Delta at centre and corners of grid cells155 148 REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2 156 149 REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12 … … 162 155 163 156 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 164 CALL wrk_alloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw, v_oce1 , v_oce2, v_ice1 )165 CALL wrk_alloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)166 CALL wrk_alloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )157 CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 ) 158 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 159 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 167 160 168 161 #if defined key_lim2 && ! defined key_lim2_vp … … 181 174 ! 182 175 !------------------------------------------------------------------------------! 183 ! 1) Ice -Snow mass (zc1), icestrength (zpresh) !176 ! 1) Ice strength (zpresh) ! 184 177 !------------------------------------------------------------------------------! 185 178 ! 186 179 ! Put every vector to 0 187 zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp 180 delta_i(:,:) = 0._wp ; 181 zpresh (:,:) = 0._wp ; 188 182 zpreshc(:,:) = 0._wp 189 183 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 190 zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 184 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 185 shear_i(:,:) = 0._wp 191 186 192 187 #if defined key_lim3 … … 198 193 !CDIR NOVERRCHK 199 194 DO ji = 1 , jpi 200 zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )201 195 #if defined key_lim3 202 196 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) … … 220 214 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 221 215 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 222 zusw(ji,jj) = 1.0 / MAX( zstms, epsd )223 216 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 224 217 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 225 218 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 226 219 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 227 & ) * zusw(ji,jj)220 & ) / MAX( zstms, epsd ) 228 221 END DO 229 222 END DO … … 267 260 DO ji = fs_2, fs_jpim1 268 261 262 zc1 = tms(ji ,jj ) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 263 zc2 = tms(ji+1,jj ) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 264 zc3 = tms(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 265 269 266 zt11 = tms(ji ,jj) * e1t(ji ,jj) 270 267 zt12 = tms(ji+1,jj) * e1t(ji+1,jj) … … 277 274 278 275 ! Mass, coriolis coeff. and currents 279 zmass1(ji,jj) = ( zt12*zc1 (ji,jj) + zt11*zc1(ji+1,jj)) / (zt11+zt12+epsd)280 zmass2(ji,jj) = ( zt22*zc1 (ji,jj) + zt21*zc1(ji,jj+1)) / (zt21+zt22+epsd)276 zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd) 277 zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd) 281 278 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) & 282 279 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd ) … … 346 343 ! 347 344 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 348 !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells345 !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 349 346 !- zds(:,:): shear on northeast corner of grid cells 350 347 ! … … 355 352 ! bugs (Martin, for Miguel). 356 353 ! 357 !- ALSO: arrays zd d, zdt, zds and delta could354 !- ALSO: arrays zdt, zds and delta could 358 355 ! be removed in the future to minimise memory demand. 359 356 ! … … 363 360 ! 364 361 ! 365 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &366 & -e2u(ji-1,jj)*u_ice(ji-1,jj) &367 & +e1v(ji,jj)*v_ice(ji,jj) &368 & -e1v(ji,jj-1)*v_ice(ji,jj-1) &369 & ) &370 & / area(ji,jj)362 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 363 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 364 & +e1v(ji,jj)*v_ice(ji,jj) & 365 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 366 & ) & 367 & / area(ji,jj) 371 368 372 369 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & … … 410 407 411 408 !- Calculate Delta at centre of grid cells 412 z zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) &409 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 413 410 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 414 411 & + e1v(ji, jj ) * u_ice2(ji,jj ) & … … 417 414 & / area(ji,jj) 418 415 419 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 420 ! MV rewriting 421 ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 422 !!gm faster to replace the line above with simply: 423 !! deltat(ji,jj) = MAX( delta, creepl ) 424 !!gm end 425 deltat(ji,jj) = delta + creepl 426 ! END MV 416 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 417 delta_i(ji,jj) = delta + creepl 427 418 !-Calculate stress tensor components zs1 and zs2 428 419 !-at centre of grid cells (see section 3.5 of CICE user's guide). 429 !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) + & 430 ! & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) & 431 ! & / ( 1._wp + alphaevp * dtotel ) 432 433 !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) - & 434 ! zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) & 435 ! & / ( 1._wp + alphaevp * ecc2 * dtotel ) 436 437 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 438 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) ) & 420 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 439 421 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 440 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta t(ji,jj) * zpresh(ji,jj) ) ) &422 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 441 423 & / ( 1._wp + dtotel ) 442 424 … … 470 452 & / ( e1f(ji,jj) * e2f(ji,jj) ) 471 453 472 deltac(ji,jj)= SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl454 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 473 455 474 456 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 475 !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) / &476 ! & ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) &477 ! & / ( 1._wp + alphaevp * ecc2 * dtotel )478 479 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)480 457 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 481 & ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj)) * zpreshc(ji,jj) ) ) &458 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 482 459 & / ( 1.0 + dtotel ) 483 460 … … 514 491 !CDIR NOVERRCHK 515 492 DO ji = fs_2, fs_jpim1 516 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 517 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 493 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 518 494 z0 = zmass1(ji,jj)/dtevp 519 495 … … 525 501 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 526 502 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 527 za*( cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))528 zcca = z0+za *cangvg529 zccb = zcorl1(ji,jj) +za*zsang503 za*(u_oce1(ji,jj)) 504 zcca = z0+za 505 zccb = zcorl1(ji,jj) 530 506 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 531 507 … … 538 514 #endif 539 515 #if defined key_bdy 540 ! clem: change u_ice and v_ice at the boundary for each iteration541 516 CALL bdy_ice_lim_dyn( 'U' ) 542 517 #endif … … 547 522 DO ji = fs_2, fs_jpim1 548 523 549 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 550 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 524 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 551 525 z0 = zmass2(ji,jj)/dtevp 552 526 ! SB modif because ocean has no slip boundary condition … … 557 531 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 558 532 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 559 za2ct(ji,jj) + za*( cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))560 zcca = z0+za *cangvg561 zccb = zcorl2(ji,jj) +za*zsang533 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 534 zcca = z0+za 535 zccb = zcorl2(ji,jj) 562 536 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 563 537 … … 570 544 #endif 571 545 #if defined key_bdy 572 ! clem: change u_ice and v_ice at the boundary for each iteration573 546 CALL bdy_ice_lim_dyn( 'V' ) 574 547 #endif … … 579 552 !CDIR NOVERRCHK 580 553 DO ji = fs_2, fs_jpim1 581 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 582 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 554 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 583 555 z0 = zmass2(ji,jj)/dtevp 584 556 ! SB modif because ocean has no slip boundary condition … … 590 562 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 591 563 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 592 za2ct(ji,jj) + za*( cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))593 zcca = z0+za *cangvg594 zccb = zcorl2(ji,jj) +za*zsang564 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 565 zcca = z0+za 566 zccb = zcorl2(ji,jj) 595 567 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 596 568 … … 603 575 #endif 604 576 #if defined key_bdy 605 ! clem: change u_ice and v_ice at the boundary for each iteration606 577 CALL bdy_ice_lim_dyn( 'V' ) 607 578 #endif … … 611 582 !CDIR NOVERRCHK 612 583 DO ji = fs_2, fs_jpim1 613 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 614 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 584 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 615 585 z0 = zmass1(ji,jj)/dtevp 616 ! SB modif because ocean has no slip boundary condition617 ! GG Bug618 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) &619 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &620 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)621 586 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 622 587 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & … … 626 591 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 627 592 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 628 za*( cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))629 zcca = z0+za *cangvg630 zccb = zcorl1(ji,jj) +za*zsang593 za*(u_oce1(ji,jj)) 594 zcca = z0+za 595 zccb = zcorl1(ji,jj) 631 596 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 632 597 END DO ! ji … … 638 603 #endif 639 604 #if defined key_bdy 640 ! clem: change u_ice and v_ice at the boundary for each iteration641 605 CALL bdy_ice_lim_dyn( 'U' ) 642 606 #endif … … 661 625 ! 4) Prevent ice velocities when the ice is thin 662 626 !------------------------------------------------------------------------------! 663 !clem : add hminrhg in the namelist664 !665 627 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 666 628 ! ocean velocity, … … 670 632 !CDIR NOVERRCHK 671 633 DO ji = fs_2, fs_jpim1 672 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )673 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )674 634 zdummy = vt_i(ji,jj) 675 635 IF ( zdummy .LE. hminrhg ) THEN … … 687 647 #endif 688 648 #if defined key_bdy 689 ! clem: change u_ice and v_ice at the boundary690 649 CALL bdy_ice_lim_dyn( 'U' ) 691 650 CALL bdy_ice_lim_dyn( 'V' ) … … 694 653 DO jj = k_j1+1, k_jpj-1 695 654 DO ji = fs_2, fs_jpim1 696 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )697 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )698 655 zdummy = vt_i(ji,jj) 699 656 IF ( zdummy .LE. hminrhg ) THEN … … 717 674 !CDIR NOVERRCHK 718 675 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi 719 !- zdd(:,:), zdt(:,:): divergence and tension at centre676 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 720 677 !- zds(:,:): shear on northeast corner of grid cells 721 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )722 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 )723 678 zdummy = vt_i(ji,jj) 724 679 IF ( zdummy .LE. hminrhg ) THEN 725 680 726 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &727 & -e2u(ji-1,jj)*u_ice(ji-1,jj) &728 & +e1v(ji,jj)*v_ice(ji,jj) &729 & -e1v(ji,jj-1)*v_ice(ji,jj-1) &730 & ) &731 & / area(ji,jj)681 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 682 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 683 & +e1v(ji,jj)*v_ice(ji,jj) & 684 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 685 & ) & 686 & / area(ji,jj) 732 687 733 688 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & … … 751 706 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 752 707 753 zdst (ji,jj)= ( e2u( ji , jj ) * v_ice1(ji ,jj ) &708 zdst = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 754 709 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 755 710 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 756 711 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 757 712 758 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 759 ! & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 760 ! & ) + creepl 761 ! MV rewriting 762 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 763 deltat(ji,jj) = delta + creepl 764 ! END MV 713 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 714 delta_i(ji,jj) = delta + creepl 765 715 766 716 ENDIF ! zdummy … … 777 727 DO jj = k_j1+1, k_jpj-1 778 728 DO ji = fs_2, fs_jpim1 779 divu_i (ji,jj) = zdd (ji,jj)780 delta_i(ji,jj) = deltat(ji,jj)781 729 ! begin TECLIM change 782 zdst (ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) &730 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 783 731 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 784 732 & + e1v( ji , jj ) * u_ice2(ji,jj) & 785 733 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 786 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst (ji,jj) * zdst(ji,jj))734 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 787 735 ! end TECLIM change 788 736 END DO … … 838 786 ! 839 787 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 840 CALL wrk_dealloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw, v_oce1 , v_oce2, v_ice1 )841 CALL wrk_dealloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)842 CALL wrk_dealloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )788 CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 ) 789 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 790 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 843 791 844 792 END SUBROUTINE lim_rhg
Note: See TracChangeset
for help on using the changeset viewer.
