Changeset 921 for trunk/NEMO/LIM_SRC_3/limrhg.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limrhg.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limrhg.F90
r888 r921 139 139 zc1 , & !: ice mass 140 140 zusw , & !: temporary weight for the computation 141 141 !: of ice strength 142 142 u_oce1, v_oce1, & !: ocean u/v component on U points 143 143 u_oce2, v_oce2, & !: ocean u/v component on V points … … 180 180 zresr !: Local error on velocity 181 181 182 !183 !------------------------------------------------------------------------------!184 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) !185 !------------------------------------------------------------------------------!186 !182 ! 183 !------------------------------------------------------------------------------! 184 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! 185 !------------------------------------------------------------------------------! 186 ! 187 187 ! Put every vector to 0 188 188 zpresh(:,:) = 0.0 ; zc1(:,:) = 0.0 … … 203 203 ! tmi = 1 where there is ice or on land 204 204 tmi(ji,jj) = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 205 205 epsd ) ) ) * tms(ji,jj) 206 206 END DO 207 207 END DO … … 213 213 !CDIR NOVERRCHK 214 214 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 215 216 217 218 219 220 221 222 223 224 215 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 216 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 217 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 218 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 219 zusw(ji,jj) = 1.0 / MAX( zstms, epsd ) 220 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 221 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 222 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 223 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 224 & ) * zusw(ji,jj) 225 225 END DO 226 226 END DO 227 227 228 228 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 229 !230 !------------------------------------------------------------------------------!231 ! 2) Wind / ocean stress, mass terms, coriolis terms232 !------------------------------------------------------------------------------!233 !229 ! 230 !------------------------------------------------------------------------------! 231 ! 2) Wind / ocean stress, mass terms, coriolis terms 232 !------------------------------------------------------------------------------! 233 ! 234 234 ! Wind stress, coriolis and mass terms on the sides of the squares 235 235 ! zfrld1: lead fraction on U-points … … 244 244 ! u_oce2: ocean u component on v points 245 245 ! v_oce2: ocean v component on v points 246 246 247 247 DO jj = k_j1+1, k_jpj-1 248 248 DO ji = fs_2, fs_jpim1 … … 255 255 ! Leads area. 256 256 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + & 257 & zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )257 & zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd ) 258 258 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + & 259 & zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )259 & zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd ) 260 260 261 261 ! Mass, coriolis coeff. and currents … … 263 263 zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) 264 264 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + & 265 266 265 e1t(ji,jj)*fcor(ji+1,jj) ) & 266 / (e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 267 267 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + & 268 269 268 e2t(ji,jj)*fcor(ji,jj+1) ) & 269 / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 270 270 ! 271 271 u_oce1(ji,jj) = u_oce(ji,jj) … … 274 274 ! Ocean has no slip boundary condition 275 275 v_oce1(ji,jj) = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj) & 276 277 276 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) & 277 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 278 278 279 279 u_oce2(ji,jj) = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj) & 280 281 280 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) & 281 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 282 282 283 283 ! Wind stress. … … 302 302 END DO 303 303 304 !305 !------------------------------------------------------------------------------!306 ! 3) Solution of the momentum equation, iterative procedure307 !------------------------------------------------------------------------------!308 !304 ! 305 !------------------------------------------------------------------------------! 306 ! 3) Solution of the momentum equation, iterative procedure 307 !------------------------------------------------------------------------------! 308 ! 309 309 ! Time step for subcycling 310 310 dtevp = rdt_ice / nevp … … 319 319 zs12(:,:) = stress12_i(:,:) 320 320 321 321 !----------------------! 322 322 DO jter = 1 , nevp ! loop over jter ! 323 323 !----------------------! 324 324 DO jj = k_j1, k_jpj-1 325 325 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 326 326 zv_ice(:,jj) = v_ice(:,jj) 327 END DO 327 END DO 328 328 329 329 DO jj = k_j1+1, k_jpj-1 330 330 DO ji = fs_2, fs_jpim1 331 331 332 !333 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)334 !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells335 !- zds(:,:): shear on northeast corner of grid cells336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 !CDIR NOVERRCHK 394 395 !CDIR NOVERRCHK 396 397 398 399 400 401 402 403 404 405 406 407 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 !CDIR NOVERRCHK 431 432 !CDIR NOVERRCHK 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 332 ! 333 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 334 !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells 335 !- zds(:,:): shear on northeast corner of grid cells 336 ! 337 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 338 ! there are many repeated calculations. 339 ! Speed could be improved by regrouping terms. For 340 ! the moment, however, the stress is on clarity of coding to avoid 341 ! bugs (Martin, for Miguel). 342 ! 343 !- ALSO: arrays zdd, zdt, zds and delta could 344 ! be removed in the future to minimise memory demand. 345 ! 346 !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 347 ! grid cells, exactly as in the B grid case. For simplicity, the indexation on 348 ! the corners is the same as in the B grid. 349 ! 350 ! 351 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 352 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 353 & +e1v(ji,jj)*v_ice(ji,jj) & 354 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 355 & ) & 356 & / area(ji,jj) 357 358 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 359 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 360 & )*e2t(ji,jj)*e2t(ji,jj) & 361 & -( v_ice(ji,jj)/e1v(ji,jj) & 362 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 363 & )*e1t(ji,jj)*e1t(ji,jj) & 364 & ) & 365 & / area(ji,jj) 366 367 ! 368 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 369 & -u_ice(ji,jj)/e1u(ji,jj) & 370 & )*e1f(ji,jj)*e1f(ji,jj) & 371 & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & 372 & -v_ice(ji,jj)/e2v(ji,jj) & 373 & )*e2f(ji,jj)*e2f(ji,jj) & 374 & ) & 375 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 376 & * tmi(ji,jj) * tmi(ji,jj+1) & 377 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 378 379 380 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 381 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 382 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 383 384 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 385 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 386 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 387 388 END DO 389 END DO 390 CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 391 CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 392 393 !CDIR NOVERRCHK 394 DO jj = k_j1+1, k_jpj-1 395 !CDIR NOVERRCHK 396 DO ji = fs_2, fs_jpim1 397 398 !- Calculate Delta at centre of grid cells 399 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & 400 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 401 & + e1v( ji , jj ) * u_ice2(ji,jj) & 402 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) & 403 & ) & 404 & / area(ji,jj) 405 406 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 407 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 408 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + & 409 (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) 410 411 !-Calculate stress tensor components zs1 and zs2 412 !-at centre of grid cells (see section 3.5 of CICE user's guide). 413 zs1(ji,jj) = ( zs1(ji,jj) & 414 & - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) + & 415 & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) & 416 * zpresh(ji,jj) ) ) & 417 & / ( 1.0 + alphaevp * dtotel ) 418 419 zs2(ji,jj) = ( zs2(ji,jj) & 420 & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) - & 421 zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) & 422 & / ( 1.0 + alphaevp*ecc2*dtotel ) 423 424 END DO 425 END DO 426 427 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 428 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 429 430 !CDIR NOVERRCHK 431 DO jj = k_j1+1, k_jpj-1 432 !CDIR NOVERRCHK 433 DO ji = fs_2, fs_jpim1 434 !- Calculate Delta on corners 435 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 436 & -v_ice1(ji,jj)/e1u(ji,jj) & 437 & )*e1f(ji,jj)*e1f(ji,jj) & 438 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 439 & -u_ice2(ji,jj)/e2v(ji,jj) & 440 & )*e2f(ji,jj)*e2f(ji,jj) & 441 & ) & 442 & / ( e1f(ji,jj) * e2f(ji,jj) ) 443 444 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 445 & -v_ice1(ji,jj)/e1u(ji,jj) & 446 & )*e1f(ji,jj)*e1f(ji,jj) & 447 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 448 & -u_ice2(ji,jj)/e2v(ji,jj) & 449 & )*e2f(ji,jj)*e2f(ji,jj) & 450 & ) & 451 & / ( e1f(ji,jj) * e2f(ji,jj) ) 452 453 deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 454 455 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 456 zs12(ji,jj) = ( zs12(ji,jj) & 457 & - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / & 458 & ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) & 459 & / ( 1.0 + alphaevp*ecc2*dtotel ) 460 461 END DO ! ji 462 END DO ! jj 463 464 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 465 466 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 467 DO jj = k_j1+1, k_jpj-1 468 DO ji = fs_2, fs_jpim1 469 !- contribution of zs1, zs2 and zs12 to zf1 470 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 471 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 472 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 473 & ) / ( e1u(ji,jj)*e2u(ji,jj) ) 474 ! contribution of zs1, zs2 and zs12 to zf2 475 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 476 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 477 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 478 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 479 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 480 END DO 481 END DO 482 482 ! 483 483 ! Computation of ice velocity … … 485 485 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 486 486 ! 487 488 489 !CDIR NOVERRCHK 490 491 !CDIR NOVERRCHK 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 !CDIR NOVERRCHK 515 516 !CDIR NOVERRCHK 517 518 519 520 521 522 523 524 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) &525 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)526 527 528 529 530 531 532 533 534 535 536 537 487 IF (MOD(jter,2).eq.0) THEN 488 489 !CDIR NOVERRCHK 490 DO jj = k_j1+1, k_jpj-1 491 !CDIR NOVERRCHK 492 DO ji = fs_2, fs_jpim1 493 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 494 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 495 z0 = zmass1(ji,jj)/dtevp 496 497 ! SB modif because ocean has no slip boundary condition 498 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 499 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 500 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 501 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 502 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 503 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 504 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 505 zcca = z0+za*cangvg 506 zccb = zcorl1(ji,jj)+za*zsang 507 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 508 509 END DO 510 END DO 511 512 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 513 514 !CDIR NOVERRCHK 515 DO jj = k_j1+1, k_jpj-1 516 !CDIR NOVERRCHK 517 DO ji = fs_2, fs_jpim1 518 519 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 520 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 521 z0 = zmass2(ji,jj)/dtevp 522 ! SB modif because ocean has no slip boundary condition 523 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 524 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 525 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 526 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 527 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 528 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 529 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 530 zcca = z0+za*cangvg 531 zccb = zcorl2(ji,jj)+za*zsang 532 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 533 534 END DO 535 END DO 536 537 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 538 538 539 539 ELSE 540 540 !CDIR NOVERRCHK 541 542 !CDIR NOVERRCHK 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 !CDIR NOVERRCHK 566 567 !CDIR NOVERRCHK 568 569 570 571 572 573 ! GG Bug574 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) &575 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &576 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 ENDIF594 595 IF(ln_ctl) THEN596 !--- Convergence test.597 DO jj = k_j1+1 , k_jpj-1598 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , &599 600 END DO601 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )602 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain603 ENDIF604 605 ! ! ==================== !541 DO jj = k_j1+1, k_jpj-1 542 !CDIR NOVERRCHK 543 DO ji = fs_2, fs_jpim1 544 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj) 545 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 546 z0 = zmass2(ji,jj)/dtevp 547 ! SB modif because ocean has no slip boundary condition 548 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 549 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 550 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 551 552 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 553 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 554 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 555 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 556 zcca = z0+za*cangvg 557 zccb = zcorl2(ji,jj)+za*zsang 558 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 559 560 END DO 561 END DO 562 563 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 564 565 !CDIR NOVERRCHK 566 DO jj = k_j1+1, k_jpj-1 567 !CDIR NOVERRCHK 568 DO ji = fs_2, fs_jpim1 569 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 570 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 571 z0 = zmass1(ji,jj)/dtevp 572 ! SB modif because ocean has no slip boundary condition 573 ! GG Bug 574 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 575 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 576 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 577 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 578 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 579 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 580 581 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 582 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 583 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 584 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 585 zcca = z0+za*cangvg 586 zccb = zcorl1(ji,jj)+za*zsang 587 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 588 END DO ! ji 589 END DO ! jj 590 591 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 592 593 ENDIF 594 595 IF(ln_ctl) THEN 596 !--- Convergence test. 597 DO jj = k_j1+1 , k_jpj-1 598 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 599 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 600 END DO 601 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 602 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 603 ENDIF 604 605 ! ! ==================== ! 606 606 END DO ! end loop over jter ! 607 607 ! ! ==================== ! 608 608 609 !610 !------------------------------------------------------------------------------!611 ! 4) Prevent ice velocities when the ice is thin612 !------------------------------------------------------------------------------!613 !609 ! 610 !------------------------------------------------------------------------------! 611 ! 4) Prevent ice velocities when the ice is thin 612 !------------------------------------------------------------------------------! 613 ! 614 614 ! If the ice thickness is below 1cm then ice velocity should equal the 615 615 ! ocean velocity, … … 636 636 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 637 637 IF ( zdummy .LE. 5.0e-2 ) THEN 638 639 640 641 642 643 644 645 638 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 639 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 640 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 641 642 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 643 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 644 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 645 ENDIF ! zdummy 646 646 END DO 647 647 END DO … … 662 662 IF ( zdummy .LE. 5.0e-2 ) THEN 663 663 664 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) &665 & -e2u(ji-1,jj)*u_ice(ji-1,jj) &666 & +e1v(ji,jj)*v_ice(ji,jj) &667 & -e1v(ji,jj-1)*v_ice(ji,jj-1) &668 & ) &669 & / area(ji,jj)670 671 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) &672 & -u_ice(ji-1,jj)/e2u(ji-1,jj) &673 & )*e2t(ji,jj)*e2t(ji,jj) &674 & -( v_ice(ji,jj)/e1v(ji,jj) &675 & -v_ice(ji,jj-1)/e1v(ji,jj-1) &676 & )*e1t(ji,jj)*e1t(ji,jj) &677 & ) &678 & / area(ji,jj)679 !680 ! SB modif because ocean has no slip boundary condition681 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) &682 & - u_ice(ji,jj) / e1u(ji,jj) ) &683 & * e1f(ji,jj) * e1f(ji,jj) &684 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) &685 & - v_ice(ji,jj) / e2v(ji,jj) ) &686 & * e2f(ji,jj) * e2f(ji,jj) ) &687 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &688 & * tmi(ji,jj) * tmi(ji,jj+1) &689 & * tmi(ji+1,jj) * tmi(ji+1,jj+1)690 691 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) &692 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) &693 & + e1v( ji , jj ) * u_ice2(ji,jj) &694 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) &695 & ) &696 & / area(ji,jj)697 698 deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + &699 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 &700 & ) + creepl701 702 664 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 665 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 666 & +e1v(ji,jj)*v_ice(ji,jj) & 667 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 668 & ) & 669 & / area(ji,jj) 670 671 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 672 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 673 & )*e2t(ji,jj)*e2t(ji,jj) & 674 & -( v_ice(ji,jj)/e1v(ji,jj) & 675 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 676 & )*e1t(ji,jj)*e1t(ji,jj) & 677 & ) & 678 & / area(ji,jj) 679 ! 680 ! SB modif because ocean has no slip boundary condition 681 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & 682 & - u_ice(ji,jj) / e1u(ji,jj) ) & 683 & * e1f(ji,jj) * e1f(ji,jj) & 684 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & 685 & - v_ice(ji,jj) / e2v(ji,jj) ) & 686 & * e2f(ji,jj) * e2f(ji,jj) ) & 687 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 688 & * tmi(ji,jj) * tmi(ji,jj+1) & 689 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 690 691 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj) & 692 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 693 & + e1v( ji , jj ) * u_ice2(ji,jj) & 694 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) & 695 & ) & 696 & / area(ji,jj) 697 698 deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 699 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & 700 & ) + creepl 701 702 ENDIF ! zdummy 703 703 704 704 END DO !jj 705 705 END DO !ji 706 !707 !------------------------------------------------------------------------------!708 ! 5) Store stress tensor and its invariants709 !------------------------------------------------------------------------------!710 !706 ! 707 !------------------------------------------------------------------------------! 708 ! 5) Store stress tensor and its invariants 709 !------------------------------------------------------------------------------! 710 ! 711 711 ! * Invariants of the stress tensor are required for limitd_me 712 712 ! accelerates convergence and improves stability … … 729 729 stress12_i(:,:) = zs12(:,:) 730 730 731 !732 !------------------------------------------------------------------------------!733 ! 6) Control prints of residual and charge ellipse734 !------------------------------------------------------------------------------!735 !731 ! 732 !------------------------------------------------------------------------------! 733 ! 6) Control prints of residual and charge ellipse 734 !------------------------------------------------------------------------------! 735 ! 736 736 ! print the residual for convergence 737 737 IF(ln_ctl) THEN
Note: See TracChangeset
for help on using the changeset viewer.
