- Timestamp:
- 2015-02-03T18:11:02+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5051 r5053 152 152 153 153 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 154 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 154 155 !!------------------------------------------------------------------- 155 156 … … 161 162 #if defined key_lim2 && ! defined key_lim2_vp 162 163 # if defined key_agrif 163 USE ice_2, vt_s => hsnm164 USE ice_2, vt_i => hicm164 USE ice_2, vt_s => hsnm 165 USE ice_2, vt_i => hicm 165 166 # else 166 vt_s => hsnm167 vt_i => hicm167 vt_s => hsnm 168 vt_i => hicm 168 169 # endif 169 at_i(:,:) = 1. - frld(:,:)170 at_i(:,:) = 1. - frld(:,:) 170 171 #endif 171 172 #if defined key_agrif && defined key_lim2 172 CALL agrif_rhg_lim2_load ! First interpolation of coarse values173 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 173 174 #endif 174 175 ! … … 189 190 #endif 190 191 191 !CDIR NOVERRCHK192 192 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 193 !CDIR NOVERRCHK194 193 DO ji = 1 , jpi 195 194 #if defined key_lim3 … … 206 205 ! Ice strength on grid cell corners (zpreshc) 207 206 ! needed for calculation of shear stress 208 !CDIR NOVERRCHK209 207 DO jj = k_j1+1, k_jpj-1 210 !CDIR NOVERRCHK211 208 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 212 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 213 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 214 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 215 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 216 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 217 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 218 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 219 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 209 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 210 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + tms(ji,jj) * wght(ji+1,jj+1,1,1) 211 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 212 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 220 213 & ) / MAX( zstms, zepsi ) 221 214 END DO 222 215 END DO 223 224 216 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 225 217 ! … … 242 234 243 235 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! 244 245 246 236 ! 237 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 238 ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 247 239 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 248 249 250 240 ! 241 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 242 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 251 243 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 252 244 ! 253 245 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 254 246 ! 255 247 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! 256 248 zpice(:,:) = ssh_m(:,:) … … 363 355 ! 364 356 ! 365 divu_i(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) 371 372 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 373 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 374 & )*e2t(ji,jj)*e2t(ji,jj) & 375 & -( v_ice(ji,jj)/e1v(ji,jj) & 376 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 377 & )*e1t(ji,jj)*e1t(ji,jj) & 378 & ) & 379 & / area(ji,jj) 357 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 358 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 359 & ) / area(ji,jj) 360 361 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 362 & - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 363 & ) / area(ji,jj) 380 364 381 365 ! 382 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 383 & -u_ice(ji,jj)/e1u(ji,jj) & 384 & )*e1f(ji,jj)*e1f(ji,jj) & 385 & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & 386 & -v_ice(ji,jj)/e2v(ji,jj) & 387 & )*e2f(ji,jj)*e2f(ji,jj) & 388 & ) & 389 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 390 & * tmi(ji,jj) * tmi(ji,jj+1) & 391 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 392 393 394 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 395 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 396 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 397 398 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 399 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 400 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 366 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 367 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 368 & ) / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2._wp - tmf(ji,jj) ) & 369 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 370 371 372 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji,jj) + v_ice(ji,jj-1) ) * e1t(ji+1,jj) & 373 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji,jj) ) & 374 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 375 376 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj) + u_ice(ji-1,jj) ) * e2t(ji,jj+1) & 377 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj) ) & 378 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 401 379 402 380 END DO … … 404 382 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 405 383 406 !CDIR NOVERRCHK407 384 DO jj = k_j1+1, k_jpj-1 408 !CDIR NOVERRCHK409 385 DO ji = fs_2, fs_jpim1 410 386 411 387 !- Calculate Delta at centre of grid cells 412 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 413 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 414 & + e1v(ji, jj ) * u_ice2(ji,jj ) & 415 & - e1v(ji, jj-1) * u_ice2(ji,jj-1) & 416 & ) & 417 & / area(ji,jj) 418 419 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 388 zdst = ( e2u(ji, jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 389 & + e1v(ji, jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 390 & ) / area(ji,jj) 391 392 delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) * usecc2 ) 420 393 delta_i(ji,jj) = delta + creepl 421 394 !-Calculate stress tensor components zs1 and zs2 … … 432 405 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 433 406 434 !CDIR NOVERRCHK435 407 DO jj = k_j1+1, k_jpj-1 436 !CDIR NOVERRCHK437 408 DO ji = fs_2, fs_jpim1 438 409 !- Calculate Delta on corners … … 490 461 IF (MOD(jter,2).eq.0) THEN 491 462 492 !CDIR NOVERRCHK493 463 DO jj = k_j1+1, k_jpj-1 494 !CDIR NOVERRCHK495 464 DO ji = fs_2, fs_jpim1 496 465 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) … … 520 489 #endif 521 490 522 !CDIR NOVERRCHK523 491 DO jj = k_j1+1, k_jpj-1 524 !CDIR NOVERRCHK525 492 DO ji = fs_2, fs_jpim1 526 493 … … 551 518 552 519 ELSE 553 !CDIR NOVERRCHK554 520 DO jj = k_j1+1, k_jpj-1 555 !CDIR NOVERRCHK556 521 DO ji = fs_2, fs_jpim1 557 522 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) … … 581 546 #endif 582 547 583 !CDIR NOVERRCHK584 548 DO jj = k_j1+1, k_jpj-1 585 !CDIR NOVERRCHK586 549 DO ji = fs_2, fs_jpim1 587 550 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) … … 628 591 ! 4) Prevent ice velocities when the ice is thin 629 592 !------------------------------------------------------------------------------! 630 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 631 ! ocean velocity, 632 ! This prevents high velocity when ice is thin 633 !CDIR NOVERRCHK 593 ! If the ice volume is below zvmin then ice velocity should equal the 594 ! ocean velocity. This prevents high velocity when ice is thin 634 595 DO jj = k_j1+1, k_jpj-1 635 !CDIR NOVERRCHK636 596 DO ji = fs_2, fs_jpim1 637 IF ( vt_i(ji,jj) .LE. hminrhg) THEN597 IF ( vt_i(ji,jj) <= zvmin ) THEN 638 598 u_ice(ji,jj) = u_oce(ji,jj) 639 599 v_ice(ji,jj) = v_oce(ji,jj) … … 655 615 DO jj = k_j1+1, k_jpj-1 656 616 DO ji = fs_2, fs_jpim1 657 IF ( vt_i(ji,jj) .LE. hminrhg) THEN658 v_ice1(ji,jj) = 0.5 *( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)&659 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))&660 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)661 662 u_ice2(ji,jj) = 0.5 *( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)&663 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))&664 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)617 IF ( vt_i(ji,jj) <= zvmin ) THEN 618 v_ice1(ji,jj) = 0.5_wp * ( ( 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 622 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 623 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 624 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 665 625 ENDIF 666 626 END DO … … 671 631 672 632 ! Recompute delta, shear and div, inputs for mechanical redistribution 673 !CDIR NOVERRCHK674 633 DO jj = k_j1+1, k_jpj-1 675 !CDIR NOVERRCHK676 634 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi 677 635 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 678 636 !- zds(:,:): shear on northeast corner of grid cells 679 IF ( vt_i(ji,jj) .LE. hminrhg) THEN637 IF ( vt_i(ji,jj) <= zvmin ) THEN 680 638 681 639 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & … … 722 680 ! 5) Store stress tensor and its invariants 723 681 !------------------------------------------------------------------------------! 724 !725 682 ! * Invariants of the stress tensor are required for limitd_me 726 683 ! (accelerates convergence and improves stability) 727 684 DO jj = k_j1+1, k_jpj-1 728 685 DO ji = fs_2, fs_jpim1 729 ! begin TECLIM change730 686 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 731 687 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & … … 733 689 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 734 690 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 735 ! end TECLIM change736 691 END DO 737 692 END DO
Note: See TracChangeset
for help on using the changeset viewer.
