- Timestamp:
- 2017-09-08T18:19:17+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8512 r8515 52 52 CONTAINS 53 53 54 SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i,delta_i )54 SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i ) 55 55 !!------------------------------------------------------------------- 56 56 !! *** SUBROUTINE ice_rhg_evp *** … … 108 108 !!------------------------------------------------------------------- 109 109 INTEGER, INTENT(in) :: kt ! time step 110 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: stress1_i, stress2_i,stress12_i111 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: u_ice, v_ice, shear_i, divu_i,delta_i110 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i 111 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 112 112 !! 113 113 INTEGER :: ji, jj ! dummy loop indices … … 247 247 248 248 ! Initialise stress tensor 249 zs1 (:,:) = stress1_i (:,:)250 zs2 (:,:) = stress2_i (:,:)251 zs12(:,:) = stress12_i(:,:)249 zs1 (:,:) = pstress1_i (:,:) 250 zs2 (:,:) = pstress2_i (:,:) 251 zs12(:,:) = pstress12_i(:,:) 252 252 253 253 ! Ice strength … … 340 340 IF(ln_ctl) THEN ! Convergence test 341 341 DO jj = 1, jpjm1 342 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step343 zv_ice(:,jj) = v_ice(:,jj)342 zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step 343 zv_ice(:,jj) = pv_ice(:,jj) 344 344 END DO 345 345 ENDIF … … 350 350 351 351 ! shear at F points 352 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) -u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &353 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) -v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &352 zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 353 & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 354 354 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 355 355 … … 367 367 368 368 ! divergence at T points 369 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) *u_ice(ji-1,jj) &370 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) *v_ice(ji,jj-1) &369 zdiv = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & 370 & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & 371 371 & ) * r1_e1e2t(ji,jj) 372 372 zdiv2 = zdiv * zdiv 373 373 374 374 ! tension at T points 375 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) -u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &376 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) -v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &375 zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 376 & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 377 377 & ) * r1_e1e2t(ji,jj) 378 378 zdt2 = zdt * zdt … … 426 426 ! 427 427 ! !--- u_ice at V point 428 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) +u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &429 & + ( u_ice(ji,jj+1) +u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)428 u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj ) + pu_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 429 & + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 430 430 ! 431 431 ! !--- v_ice at U point 432 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) +v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &433 & + ( v_ice(ji+1,jj) +v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)432 v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji ,jj) + pv_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 433 & + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 434 434 ! 435 435 END DO … … 444 444 DO ji = fs_2, fs_jpim1 445 445 ! !--- tau_io/(v_oce - v_ice) 446 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &446 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & 447 447 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 448 448 ! !--- Ocean-to-Ice stress 449 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )449 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 450 450 ! 451 451 ! !--- tau_bottom/v_ice 452 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) *v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )452 zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 453 453 zTauB = - tau_icebfr(ji,jj) / zvel 454 454 ! 455 455 ! !--- Coriolis at V-points (energy conserving formulation) 456 456 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 457 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) *u_ice(ji-1,jj ) ) &458 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) *u_ice(ji-1,jj+1) ) )457 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & 458 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 459 459 ! 460 460 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 465 465 ! 466 466 ! !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 467 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj)& ! previous velocity468 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)469 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast470 & + ( 1._wp - rswitch ) *v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0471 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )& ! v_ice = v_oce if mass < zmmin472 & ) * zmaskV(ji,jj)467 pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity 468 & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 469 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 470 & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 471 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 472 & ) * zmaskV(ji,jj) 473 473 ! 474 474 ! Bouillon 2013 475 !! v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) +v_ice_b(ji,jj) ) &475 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & 476 476 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 477 477 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) … … 479 479 END DO 480 480 END DO 481 CALL lbc_lnk( v_ice, 'V', -1. )481 CALL lbc_lnk( pv_ice, 'V', -1. ) 482 482 ! 483 483 #if defined key_agrif … … 491 491 492 492 ! tau_io/(u_oce - u_ice) 493 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &493 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & 494 494 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 495 495 496 496 ! Ocean-to-Ice stress 497 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )497 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 498 498 499 499 ! tau_bottom/u_ice 500 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) *u_ice(ji,jj) ) )500 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 501 501 zTauB = - tau_icebfr(ji,jj) / zvel 502 502 503 503 ! Coriolis at U-points (energy conserving formulation) 504 504 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 505 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) *v_ice(ji ,jj-1) ) &506 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) *v_ice(ji+1,jj-1) ) )505 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & 506 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 507 507 508 508 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 513 513 514 514 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 515 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj)& ! previous velocity516 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)517 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast518 & + ( 1._wp - rswitch ) *u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0519 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )& ! v_ice = v_oce if mass < zmmin520 & ) * zmaskU(ji,jj)515 pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity 516 & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 517 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 518 & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 519 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 520 & ) * zmaskU(ji,jj) 521 521 ! Bouillon 2013 522 !! u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) +u_ice_b(ji,jj) ) &522 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & 523 523 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 524 524 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 525 525 END DO 526 526 END DO 527 CALL lbc_lnk( u_ice, 'U', -1. )527 CALL lbc_lnk( pu_ice, 'U', -1. ) 528 528 ! 529 529 #if defined key_agrif … … 539 539 540 540 ! tau_io/(u_oce - u_ice) 541 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &541 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & 542 542 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 543 543 544 544 ! Ocean-to-Ice stress 545 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )545 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 546 546 547 547 ! tau_bottom/u_ice 548 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) *u_ice(ji,jj) ) )548 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 549 549 zTauB = - tau_icebfr(ji,jj) / zvel 550 550 551 551 ! Coriolis at U-points (energy conserving formulation) 552 552 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 553 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) *v_ice(ji ,jj-1) ) &554 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) *v_ice(ji+1,jj-1) ) )553 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & 554 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 555 555 556 556 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 561 561 562 562 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 563 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj)& ! previous velocity564 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)565 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast566 & + ( 1._wp - rswitch ) *u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0567 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )& ! v_ice = v_oce if mass < zmmin568 & ) * zmaskU(ji,jj)563 pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity 564 & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 565 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 566 & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 567 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 568 & ) * zmaskU(ji,jj) 569 569 ! Bouillon 2013 570 !! u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) +u_ice_b(ji,jj) ) &570 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & 571 571 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 572 572 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 573 573 END DO 574 574 END DO 575 CALL lbc_lnk( u_ice, 'U', -1. )575 CALL lbc_lnk( pu_ice, 'U', -1. ) 576 576 ! 577 577 #if defined key_agrif … … 585 585 586 586 ! tau_io/(v_oce - v_ice) 587 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &587 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & 588 588 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 589 589 590 590 ! Ocean-to-Ice stress 591 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )591 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 592 592 593 593 ! tau_bottom/v_ice 594 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) *v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )594 zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 595 595 ztauB = - tau_icebfr(ji,jj) / zvel 596 596 597 597 ! Coriolis at V-points (energy conserving formulation) 598 598 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 599 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) *u_ice(ji-1,jj ) ) &600 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) *u_ice(ji-1,jj+1) ) )599 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & 600 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 601 601 602 602 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 607 607 608 608 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 609 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj)& ! previous velocity610 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)611 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast612 & + ( 1._wp - rswitch ) *v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0613 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )& ! v_ice = v_oce if mass < zmmin614 & ) * zmaskV(ji,jj)609 pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity 610 & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 611 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 612 & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 613 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 614 & ) * zmaskV(ji,jj) 615 615 ! Bouillon 2013 616 !! v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) +v_ice_b(ji,jj) ) &616 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & 617 617 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 618 618 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 619 619 END DO 620 620 END DO 621 CALL lbc_lnk( v_ice, 'V', -1. )621 CALL lbc_lnk( pv_ice, 'V', -1. ) 622 622 ! 623 623 #if defined key_agrif … … 631 631 IF(ln_ctl) THEN ! Convergence test 632 632 DO jj = 2 , jpjm1 633 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS(v_ice(:,jj) - zv_ice(:,jj) ) )633 zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) ) 634 634 END DO 635 635 zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) … … 648 648 649 649 ! shear at F points 650 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) -u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &651 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) -v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &650 zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 651 & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 652 652 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 653 653 … … 660 660 661 661 ! tension**2 at T points 662 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) -u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &663 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) -v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &662 zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 663 & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 664 664 & ) * r1_e1e2t(ji,jj) 665 665 zdt2 = zdt * zdt … … 671 671 672 672 ! shear at T points 673 shear_i(ji,jj) = SQRT( zdt2 + zds2 )673 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 674 674 675 675 ! divergence at T points 676 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) *u_ice(ji-1,jj) &677 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) *v_ice(ji,jj-1) &678 & ) * r1_e1e2t(ji,jj)676 pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & 677 & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & 678 & ) * r1_e1e2t(ji,jj) 679 679 680 680 ! delta at T points 681 zdelta = SQRT( divu_i(ji,jj) *divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )681 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 682 682 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 683 delta_i(ji,jj) = zdelta + rn_creepl * rswitch683 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 684 684 685 685 END DO 686 686 END DO 687 CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1.,delta_i, 'T', 1. )687 CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 688 688 689 689 ! --- Store the stress tensor for the next time step --- ! 690 stress1_i (:,:) = zs1 (:,:)691 stress2_i (:,:) = zs2 (:,:)692 stress12_i(:,:) = zs12(:,:)690 pstress1_i (:,:) = zs1 (:,:) 691 pstress2_i (:,:) = zs2 (:,:) 692 pstress12_i(:,:) = zs12(:,:) 693 693 ! 694 694 … … 703 703 704 704 ! --- divergence, shear and strength --- ! 705 IF( iom_use('idive') ) CALL iom_put( "idive" , divu_i(:,:) * zswi(:,:) ) ! divergence706 IF( iom_use('ishear') ) CALL iom_put( "ishear" , shear_i(:,:) * zswi(:,:) ) ! shear705 IF( iom_use('idive') ) CALL iom_put( "idive" , pdivu_i(:,:) * zswi(:,:) ) ! divergence 706 IF( iom_use('ishear') ) CALL iom_put( "ishear" , pshear_i(:,:) * zswi(:,:) ) ! shear 707 707 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Ice strength 708 708 … … 714 714 DO jj = 2, jpjm1 715 715 DO ji = 2, jpim1 716 zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji ,jj-1) *stress12_i(ji ,jj-1) + & ! stress12_i at T-point717 & zswi(ji ,jj) * stress12_i(ji ,jj) + zswi(ji-1,jj-1) *stress12_i(ji-1,jj-1) ) &718 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )719 720 zshear = SQRT( stress2_i(ji,jj) *stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress716 zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 717 & zswi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 718 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 719 720 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 721 721 722 722 zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 723 723 724 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)725 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)726 !! zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress724 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 725 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 726 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 727 727 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 728 zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015728 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 729 729 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 730 zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )730 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 731 731 END DO 732 732 END DO … … 773 773 774 774 ! 2D ice mass, snow mass, area transport arrays (X, Y) 775 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch776 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch775 zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch 776 zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch 777 777 778 778 zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
Note: See TracChangeset
for help on using the changeset viewer.
