Changeset 8678
- Timestamp:
- 2017-11-08T19:17:04+01:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Files:
-
- 4 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/icedyn_rhg.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/icedyn_rhg_evp.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r8598 r8678 90 90 !------------------------------------------------------------------------------ 91 91 ln_rhg_EVP = .true. ! EVP rheology 92 rn_creepl = 1.0e-12 ! creep limit [1/s] 92 ln_aEVP = .false. ! adaptive rheology (Kimmritz et al. 2016 & 2017) 93 rn_creepl = 2.0e-9 ! creep limit [1/s] 93 94 rn_ecc = 2.0 ! eccentricity of the elliptical yield curve 94 95 nn_nevp = 120 ! number of EVP subcycles -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8597 r8678 178 178 ! 179 179 ! !!** ice-rheology namelist (namrhg) ** 180 LOGICAL , PUBLIC :: ln_aEVP !: using adaptive EVP (T or F) 180 181 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 181 182 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90
r8534 r8678 78 78 CASE( np_rhgEVP ) ! Elasto-Viscous-Plastic ! 79 79 ! !------------------------! 80 CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice,shear_i, divu_i, delta_i )80 CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 81 81 ! 82 82 END SELECT … … 107 107 INTEGER :: ios, ioptio ! Local integer output status for namelist read 108 108 !! 109 NAMELIST/namdyn_rhg/ ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast109 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 110 110 !!------------------------------------------------------------------- 111 111 ! … … 125 125 WRITE(numout,*) ' Namelist namdyn_rhg:' 126 126 WRITE(numout,*) ' rheology EVP (icedyn_rhg_evp) ln_rhg_EVP = ', ln_rhg_EVP 127 WRITE(numout,*) ' use adaptive EVP (aEVP) ln_aEVP = ', ln_aEVP 127 128 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 128 129 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90
r8578 r8678 53 53 CONTAINS 54 54 55 SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice,pshear_i, pdivu_i, pdelta_i )55 SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) 56 56 !!------------------------------------------------------------------- 57 57 !! *** SUBROUTINE ice_dyn_rhg_evp *** … … 95 95 !! 5) Diagnostics including charge ellipse 96 96 !! 97 !! ** Notes : There is the possibility to use mEVP from Bouillon 2013 98 !! (by uncommenting some lines in part 3 and changing alpha and beta parameters) 99 !! but this solution appears very unstable (see Kimmritz et al 2016) 97 !! ** Notes : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017) 98 !! by setting up ln_aEVP=T (i.e. changing alpha and beta parameters). 99 !! This is an upgraded version of mEVP from Bouillon et al. 2013 100 !! (i.e. more stable and better convergence) 100 101 !! 101 102 !! References : Hunke and Dukowicz, JPO97 102 103 !! Bouillon et al., Ocean Modelling 2009 103 104 !! Bouillon et al., Ocean Modelling 2013 105 !! Kimmritz et al., Ocean Modelling 2016 & 2017 104 106 !!------------------------------------------------------------------- 105 107 INTEGER, INTENT(in) :: kt ! time step 106 108 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i 107 REAL(wp), DIMENSION(:,:), INTENT(inout) :: u_ice, v_ice108 109 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i, pdivu_i, pdelta_i 109 110 !! … … 114 115 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 115 116 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 116 REAL(wp) :: z beta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013117 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 117 118 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass 118 119 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars … … 126 127 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors 127 128 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 128 ! 129 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 130 ! 131 REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points 129 132 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 130 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! ice/snow mass/dton U/V points133 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 131 134 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 132 135 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points … … 231 234 232 235 ! alpha parameters (Bouillon 2009) 233 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 234 zalph2 = zalph1 * z1_ecc2 235 236 ! alpha and beta parameters (Bouillon 2013) 237 !!zalph1 = 40. 238 !!zalph2 = 40. 239 !!zbeta = 3000. 240 !!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001) 241 242 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 243 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 244 236 IF( .NOT. ln_aEVP ) THEN 237 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 238 zalph2 = zalph1 * z1_ecc2 239 240 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 241 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 242 ENDIF 243 245 244 ! Initialise stress tensor 246 245 zs1 (:,:) = pstress1_i (:,:) … … 304 303 zmf(ji,jj) = zm1 * ff_t(ji,jj) 305 304 305 ! dt/m at T points (for alpha and beta coefficients) 306 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 307 306 308 ! m/dt 307 309 zmU_t(ji,jj) = zmassU * z1_dtevp 308 310 zmV_t(ji,jj) = zmassV * z1_dtevp 309 311 310 312 ! Drag ice-atm. 311 313 zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) … … 326 328 END DO 327 329 END DO 328 CALL lbc_lnk ( zmf, 'T', 1. )330 CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 329 331 ! 330 332 !------------------------------------------------------------------------------! … … 380 382 ! P/delta at T points 381 383 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 384 385 ! alpha & beta for aEVP 386 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 387 ! alpha = beta = sqrt(4*gamma) 388 IF( ln_aEVP ) THEN 389 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 390 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 391 zalph2 = zalph1 392 z1_alph2 = z1_alph1 393 ENDIF 382 394 383 395 ! stress at T points … … 392 404 DO ji = 1, jpim1 393 405 406 ! alpha & beta for aEVP 407 IF( ln_aEVP ) THEN 408 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 409 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 410 zbeta(ji,jj) = zalph2 411 ENDIF 412 394 413 ! P/delta at F points 395 414 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) … … 411 430 & ) * 2._wp * r1_e1u(ji,jj) & 412 431 & ) * r1_e1e2u(ji,jj) 413 414 432 ! 433 ! !--- V points 415 434 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 416 435 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 419 438 & ) * 2._wp * r1_e2v(ji,jj) & 420 439 & ) * r1_e1e2v(ji,jj) 421 422 440 ! 441 ! !--- u_ice at V point 423 442 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 424 443 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 425 426 444 ! 445 ! !--- v_ice at U point 427 446 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 428 447 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 429 !430 448 END DO 431 449 END DO 432 450 ! 433 451 ! --- Computation of ice velocity --- ! 434 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp452 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017 435 453 ! Bouillon et al. 2009 (eq 34-35) => stable 436 454 IF( MOD(jter,2) == 0 ) THEN ! even iterations … … 438 456 DO jj = 2, jpjm1 439 457 DO ji = fs_2, fs_jpim1 440 458 ! !--- tau_io/(v_oce - v_ice) 441 459 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 442 460 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 443 461 ! !--- Ocean-to-Ice stress 444 462 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 445 446 463 ! 464 ! !--- tau_bottom/v_ice 447 465 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 448 466 zTauB = - tau_icebfr(ji,jj) / zvel 449 450 467 ! 468 ! !--- Coriolis at V-points (energy conserving formulation) 451 469 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 452 470 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 453 471 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 454 455 472 ! 473 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 456 474 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 457 458 475 ! 476 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 459 477 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 460 ! 461 ! !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 478 ! 479 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 480 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 481 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 482 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 483 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 484 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 485 & ) * zmaskV(ji,jj) 486 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 462 487 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 463 488 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 466 491 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 467 492 & ) * zmaskV(ji,jj) 468 ! 469 ! Bouillon 2013 470 !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 471 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 472 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 473 ! 493 ENDIF 474 494 END DO 475 495 END DO … … 483 503 ! 484 504 DO jj = 2, jpjm1 485 DO ji = fs_2, fs_jpim1 486 487 ! tau_io/(u_oce - u_ice) 505 DO ji = fs_2, fs_jpim1 506 ! !--- tau_io/(u_oce - u_ice) 488 507 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 489 508 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 490 491 ! Ocean-to-Ice stress 509 ! !--- Ocean-to-Ice stress 492 510 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 493 494 ! tau_bottom/u_ice511 ! 512 ! !--- tau_bottom/u_ice 495 513 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 496 514 zTauB = - tau_icebfr(ji,jj) / zvel 497 498 ! Coriolis at U-points (energy conserving formulation)515 ! 516 ! !--- Coriolis at U-points (energy conserving formulation) 499 517 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 500 518 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 501 519 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 502 503 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io520 ! 521 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 504 522 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 505 506 ! landfast switch => 0 = static friction ; 1 = sliding friction523 ! 524 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 507 525 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 508 509 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 526 ! 527 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 528 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 529 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 530 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 531 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 532 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 533 & ) * zmaskU(ji,jj) 534 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 510 535 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 511 536 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 514 539 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 515 540 & ) * zmaskU(ji,jj) 516 ! Bouillon 2013 517 !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 518 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 519 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 541 ENDIF 520 542 END DO 521 543 END DO … … 532 554 DO jj = 2, jpjm1 533 555 DO ji = fs_2, fs_jpim1 534 535 ! tau_io/(u_oce - u_ice) 556 ! !--- tau_io/(u_oce - u_ice) 536 557 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 537 558 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 538 539 ! Ocean-to-Ice stress 559 ! !--- Ocean-to-Ice stress 540 560 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 541 542 ! tau_bottom/u_ice561 ! 562 ! !--- tau_bottom/u_ice 543 563 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 544 564 zTauB = - tau_icebfr(ji,jj) / zvel 545 546 ! Coriolis at U-points (energy conserving formulation)565 ! 566 ! !--- Coriolis at U-points (energy conserving formulation) 547 567 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 548 568 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 549 569 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 550 551 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io570 ! 571 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 552 572 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 553 554 ! landfast switch => 0 = static friction ; 1 = sliding friction573 ! 574 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 555 575 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 556 557 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 576 ! 577 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 578 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 579 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 580 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 581 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 582 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 583 & ) * zmaskU(ji,jj) 584 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 558 585 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 559 586 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 560 587 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 561 588 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 562 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 589 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 563 590 & ) * zmaskU(ji,jj) 564 ! Bouillon 2013 565 !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 566 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 567 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 591 ENDIF 568 592 END DO 569 593 END DO … … 578 602 DO jj = 2, jpjm1 579 603 DO ji = fs_2, fs_jpim1 580 581 ! tau_io/(v_oce - v_ice) 604 ! !--- tau_io/(v_oce - v_ice) 582 605 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 583 606 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 584 585 ! Ocean-to-Ice stress 607 ! !--- Ocean-to-Ice stress 586 608 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 587 588 ! tau_bottom/v_ice609 ! 610 ! !--- tau_bottom/v_ice 589 611 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 590 612 ztauB = - tau_icebfr(ji,jj) / zvel 591 592 ! Coriolis at V-points (energy conserving formulation)613 ! 614 ! !--- Coriolis at v-points (energy conserving formulation) 593 615 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 594 616 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 595 617 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 596 597 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io618 ! 619 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 598 620 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 599 600 ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction621 ! 622 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 601 623 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 602 603 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 624 ! 625 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 626 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 627 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 628 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 629 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 630 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 631 & ) * zmaskV(ji,jj) 632 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 604 633 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 605 634 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 608 637 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 609 638 & ) * zmaskV(ji,jj) 610 ! Bouillon 2013 611 !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 612 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 613 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 639 ENDIF 614 640 END DO 615 641 END DO
Note: See TracChangeset
for help on using the changeset viewer.
