- Timestamp:
- 2017-11-24T17:56:51+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90
r8637 r8813 16 16 !! 'key_lim3' ESIM sea-ice model 17 17 !!---------------------------------------------------------------------- 18 !! ice_dyn_rhg_evp : computes ice velocities from EVP rheology 18 !! ice_dyn_rhg_evp : computes ice velocities from EVP rheology 19 !! rhg_evp_rst : read/write EVP fields in ice restart 19 20 !!---------------------------------------------------------------------- 20 21 USE phycst ! Physical constant … … 24 25 USE ice ! sea-ice: ice variables 25 26 USE icedyn_rdgrft ! sea-ice: ice strength 27 USE bdy_oce , ONLY : ln_bdy 28 USE bdyice 29 #if defined key_agrif 30 USE agrif_lim3_interp 31 #endif 26 32 ! 27 33 USE in_out_manager ! I/O manager … … 31 37 USE lbclnk ! lateral boundary conditions (or mpp links) 32 38 USE prtctl ! Print control 33 !34 #if defined key_agrif35 USE agrif_lim3_interp36 #endif37 USE bdy_oce , ONLY: ln_bdy38 USE bdyice39 39 40 40 IMPLICIT NONE … … 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 *** 58 !! EVP-C-grid58 !! EVP-C-grid 59 59 !! 60 60 !! ** purpose : determines sea ice drift from wind stress, ice-ocean … … 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 INTEGER, INTENT(in) :: kt ! time step 106 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i 107 REAL(wp), DIMENSION(:,:), INTENT(inout) :: u_ice, v_ice 108 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i, pdivu_i, pdelta_i 107 INTEGER , INTENT(in ) :: kt ! time step 108 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i ! 109 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 109 110 !! 110 111 INTEGER :: ji, jj ! dummy loop indices 111 112 INTEGER :: jter ! local integers 112 113 ! 113 114 REAL(wp) :: zrhoco ! rau0 * rn_cio 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 119 120 REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars 120 121 ! 121 122 REAL(wp) :: zresm ! Maximal error on ice velocity 122 123 REAL(wp) :: zintb, zintn ! dummy argument 123 124 REAL(wp) :: zfac_x, zfac_y 124 125 REAL(wp) :: zshear, zdum1, zdum2 125 126 ! 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 … … 134 137 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 135 138 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 136 139 ! 137 140 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 138 141 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 139 142 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence 140 143 REAL(wp), DIMENSION(jpi,jpj) :: zpice ! array used for the calculation of ice surface slope: 141 142 ! ice top surface if ice is embedded144 ! ! ocean surface (ssh_m) if ice is not embedded 145 ! ! ice top surface if ice is embedded 143 146 REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array 144 147 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array 145 148 ! 146 149 REAL(wp), DIMENSION(jpi,jpj) :: zswitchU, zswitchV ! dummy arrays 147 150 REAL(wp), DIMENSION(jpi,jpj) :: zmaskU, zmaskV ! mask for ice presence … … 179 182 #endif 180 183 ! 184 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 181 185 !------------------------------------------------------------------------------! 182 186 ! 0) mask at F points for the ice … … 231 235 232 236 ! 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 237 IF( .NOT. ln_aEVP ) THEN 238 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 239 zalph2 = zalph1 * z1_ecc2 240 241 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 242 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 243 ENDIF 244 245 245 ! Initialise stress tensor 246 246 zs1 (:,:) = pstress1_i (:,:) … … 266 266 IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==! 267 267 ! 268 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}269 ! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}268 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 269 ! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1} 270 270 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 271 271 ! 272 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) *{SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}272 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 273 273 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 274 274 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp … … 304 304 zmf(ji,jj) = zm1 * ff_t(ji,jj) 305 305 306 ! dt/m at T points (for alpha and beta coefficients) 307 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 308 306 309 ! m/dt 307 310 zmU_t(ji,jj) = zmassU * z1_dtevp 308 311 zmV_t(ji,jj) = zmassV * z1_dtevp 309 312 310 313 ! Drag ice-atm. 311 314 zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) … … 326 329 END DO 327 330 END DO 328 CALL lbc_lnk ( zmf, 'T', 1. )331 CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 329 332 ! 330 333 !------------------------------------------------------------------------------! … … 380 383 ! P/delta at T points 381 384 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 385 386 ! alpha & beta for aEVP 387 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 388 ! alpha = beta = sqrt(4*gamma) 389 IF( ln_aEVP ) THEN 390 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 391 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 392 zalph2 = zalph1 393 z1_alph2 = z1_alph1 394 ENDIF 382 395 383 396 ! stress at T points … … 392 405 DO ji = 1, jpim1 393 406 407 ! alpha & beta for aEVP 408 IF( ln_aEVP ) THEN 409 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 410 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 411 zbeta(ji,jj) = zalph2 412 ENDIF 413 394 414 ! P/delta at F points 395 415 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 431 & ) * 2._wp * r1_e1u(ji,jj) & 412 432 & ) * r1_e1e2u(ji,jj) 413 414 433 ! 434 ! !--- V points 415 435 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 416 436 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 419 439 & ) * 2._wp * r1_e2v(ji,jj) & 420 440 & ) * r1_e1e2v(ji,jj) 421 422 441 ! 442 ! !--- u_ice at V point 423 443 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 424 444 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 425 426 445 ! 446 ! !--- v_ice at U point 427 447 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 428 448 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 429 !430 449 END DO 431 450 END DO 432 451 ! 433 452 ! --- Computation of ice velocity --- ! 434 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp453 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017 435 454 ! Bouillon et al. 2009 (eq 34-35) => stable 436 455 IF( MOD(jter,2) == 0 ) THEN ! even iterations … … 438 457 DO jj = 2, jpjm1 439 458 DO ji = fs_2, fs_jpim1 440 459 ! !--- tau_io/(v_oce - v_ice) 441 460 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 442 461 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 443 462 ! !--- Ocean-to-Ice stress 444 463 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 445 446 464 ! 465 ! !--- tau_bottom/v_ice 447 466 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 448 467 zTauB = - tau_icebfr(ji,jj) / zvel 449 450 468 ! 469 ! !--- Coriolis at V-points (energy conserving formulation) 451 470 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 452 471 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 453 472 & + 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 473 ! 474 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 456 475 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 457 458 476 ! 477 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 459 478 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) 479 ! 480 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 481 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 482 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 483 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 484 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 485 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 486 & ) * zmaskV(ji,jj) 487 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 462 488 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 463 489 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 466 492 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 467 493 & ) * 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 ! 494 ENDIF 474 495 END DO 475 496 END DO … … 483 504 ! 484 505 DO jj = 2, jpjm1 485 DO ji = fs_2, fs_jpim1 486 487 ! tau_io/(u_oce - u_ice) 506 DO ji = fs_2, fs_jpim1 507 ! !--- tau_io/(u_oce - u_ice) 488 508 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 489 509 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 490 491 ! Ocean-to-Ice stress 510 ! !--- Ocean-to-Ice stress 492 511 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 493 494 ! tau_bottom/u_ice512 ! 513 ! !--- tau_bottom/u_ice 495 514 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 496 515 zTauB = - tau_icebfr(ji,jj) / zvel 497 498 ! Coriolis at U-points (energy conserving formulation)516 ! 517 ! !--- Coriolis at U-points (energy conserving formulation) 499 518 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 500 519 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 501 520 & + 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_io521 ! 522 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 504 523 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 friction524 ! 525 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 507 526 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) 527 ! 528 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 529 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 530 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 531 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 532 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 533 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 534 & ) * zmaskU(ji,jj) 535 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 510 536 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 511 537 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 514 540 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 515 541 & ) * 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) 542 ENDIF 520 543 END DO 521 544 END DO … … 532 555 DO jj = 2, jpjm1 533 556 DO ji = fs_2, fs_jpim1 534 535 ! tau_io/(u_oce - u_ice) 557 ! !--- tau_io/(u_oce - u_ice) 536 558 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 537 559 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 538 539 ! Ocean-to-Ice stress 560 ! !--- Ocean-to-Ice stress 540 561 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 541 542 ! tau_bottom/u_ice562 ! 563 ! !--- tau_bottom/u_ice 543 564 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 544 565 zTauB = - tau_icebfr(ji,jj) / zvel 545 546 ! Coriolis at U-points (energy conserving formulation)566 ! 567 ! !--- Coriolis at U-points (energy conserving formulation) 547 568 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 548 569 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 549 570 & + 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_io571 ! 572 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 552 573 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 friction574 ! 575 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 555 576 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) 577 ! 578 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 579 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 580 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 581 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 582 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 583 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 584 & ) * zmaskU(ji,jj) 585 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 558 586 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 559 587 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 560 588 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 561 589 & + ( 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 590 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 563 591 & ) * 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) 592 ENDIF 568 593 END DO 569 594 END DO … … 578 603 DO jj = 2, jpjm1 579 604 DO ji = fs_2, fs_jpim1 580 581 ! tau_io/(v_oce - v_ice) 605 ! !--- tau_io/(v_oce - v_ice) 582 606 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 583 607 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 584 585 ! Ocean-to-Ice stress 608 ! !--- Ocean-to-Ice stress 586 609 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 587 588 ! tau_bottom/v_ice610 ! 611 ! !--- tau_bottom/v_ice 589 612 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 590 613 ztauB = - tau_icebfr(ji,jj) / zvel 591 592 ! Coriolis at V-points (energy conserving formulation)614 ! 615 ! !--- Coriolis at v-points (energy conserving formulation) 593 616 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 594 617 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 595 618 & + 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_io619 ! 620 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 598 621 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 friction622 ! 623 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 601 624 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) 625 ! 626 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 627 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 628 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 629 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 630 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 631 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 632 & ) * zmaskV(ji,jj) 633 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 604 634 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 605 635 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 608 638 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 609 639 & ) * 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) 640 ENDIF 614 641 END DO 615 642 END DO … … 820 847 END SUBROUTINE ice_dyn_rhg_evp 821 848 849 822 850 SUBROUTINE rhg_evp_rst( cdrw, kt ) 823 851 !!---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
