Changeset 10413 for NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2018-12-18T18:59:59+01:00 (6 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r10332 r10413 119 119 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 120 120 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 121 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass121 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 122 122 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 123 123 REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars 124 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 125 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 124 126 ! 125 127 REAL(wp) :: zresm ! Maximal error on ice velocity … … 137 139 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 138 140 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points 141 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ib , ztauV_ib ! ice-bottom stress at U-V points (landfast param) 139 142 REAL(wp), DIMENSION(jpi,jpj) :: zspgU , zspgV ! surface pressure gradient at U/V points 140 143 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points … … 256 259 END DO 257 260 END DO 258 261 262 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 263 IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN ; zkt = rn_tensile 264 ELSE ; zkt = 0._wp 265 ENDIF 259 266 ! 260 267 !------------------------------------------------------------------------------! … … 317 324 CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 318 325 ! 326 ! !== Landfast ice parameterization ==! 327 ! 328 IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 329 DO jj = 2, jpjm1 330 DO ji = fs_2, fs_jpim1 331 ! ice thickness at U-V points 332 zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 333 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 334 ! ice-bottom stress at U points 335 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 336 zTauU_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 337 ! ice-bottom stress at V points 338 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 339 zTauV_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 340 ! ice_bottom stress at T points 341 zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 342 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 343 END DO 344 END DO 345 CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. ) 346 ! 347 ELSEIF( ln_landfast_home ) THEN !-- Home made 348 DO jj = 2, jpjm1 349 DO ji = fs_2, fs_jpim1 350 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 351 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 352 END DO 353 END DO 354 ! 355 ELSE !-- no landfast 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 358 zTauU_ib(ji,jj) = 0._wp 359 zTauV_ib(ji,jj) = 0._wp 360 END DO 361 END DO 362 ENDIF 363 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 364 319 365 !------------------------------------------------------------------------------! 320 366 ! 3) Solution of the momentum equation, iterative procedure … … 380 426 ENDIF 381 427 382 ! stress at T points 383 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta) ) * z1_alph1384 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2428 ! stress at T points (zkt/=0 if landfast) 429 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 430 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 385 431 386 432 END DO … … 401 447 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) ) 402 448 403 ! stress at F points 404 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2449 ! stress at F points (zkt/=0 if landfast) 450 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 405 451 406 452 END DO … … 450 496 ! 451 497 ! !--- 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)) )453 zTauB = - tau_icebfr(ji,jj) / zvel498 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 499 zTauB = - zTauV_ib(ji,jj) / zvel 454 500 ! 455 501 ! !--- Coriolis at V-points (energy conserving formulation) … … 462 508 ! 463 509 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 464 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )510 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 465 511 ! 466 512 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 498 544 ! 499 545 ! !--- 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)) )501 zTauB = - tau_icebfr(ji,jj) / zvel546 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 547 zTauB = - zTauU_ib(ji,jj) / zvel 502 548 ! 503 549 ! !--- Coriolis at U-points (energy conserving formulation) … … 510 556 ! 511 557 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 512 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )558 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 513 559 ! 514 560 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 548 594 ! 549 595 ! !--- tau_bottom/u_ice 550 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj)) )551 zTauB = - tau_icebfr(ji,jj) / zvel596 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 597 zTauB = - zTauU_ib(ji,jj) / zvel 552 598 ! 553 599 ! !--- Coriolis at U-points (energy conserving formulation) … … 560 606 ! 561 607 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 562 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )608 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 563 609 ! 564 610 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 596 642 ! 597 643 ! !--- tau_bottom/v_ice 598 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj)) )599 z tauB = - tau_icebfr(ji,jj) / zvel644 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 645 zTauB = - zTauV_ib(ji,jj) / zvel 600 646 ! 601 647 ! !--- Coriolis at v-points (energy conserving formulation) … … 608 654 ! 609 655 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 610 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )656 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 611 657 ! 612 658 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
Note: See TracChangeset
for help on using the changeset viewer.
