- Timestamp:
- 2018-11-15T10:46:57+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rhg_evp.F90
r9935 r10312 118 118 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 119 119 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 120 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass120 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 121 121 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 122 122 REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars 123 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 124 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 123 125 ! 124 126 REAL(wp) :: zresm ! Maximal error on ice velocity … … 136 138 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 137 139 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points 140 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ib , ztauV_ib ! ice-bottom stress at U-V points (landfast param) 138 141 REAL(wp), DIMENSION(jpi,jpj) :: zspgU , zspgV ! surface pressure gradient at U/V points 139 142 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points … … 255 258 END DO 256 259 END DO 257 260 261 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 262 IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN ; zkt = rn_tensile 263 ELSE ; zkt = 0._wp 264 ENDIF 258 265 ! 259 266 !------------------------------------------------------------------------------! … … 273 280 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 274 281 ! 275 ELSE 282 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! 276 283 zpice(:,:) = ssh_m(:,:) 277 284 ENDIF … … 328 335 CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 329 336 ! 337 ! !== Landfast ice parameterization ==! 338 ! 339 IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 340 DO jj = 2, jpjm1 341 DO ji = fs_2, fs_jpim1 342 ! ice thickness at U-V points 343 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) 344 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) 345 ! ice-bottom stress at U points 346 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 347 zTauU_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 348 ! ice-bottom stress at V points 349 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 350 zTauV_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 351 ! ice_bottom stress at T points 352 zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 353 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 354 END DO 355 END DO 356 CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. ) 357 ! 358 ELSEIF( ln_landfast_home ) THEN !-- Home made 359 DO jj = 2, jpjm1 360 DO ji = fs_2, fs_jpim1 361 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 362 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 363 END DO 364 END DO 365 ! 366 ELSE !-- no landfast 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 369 zTauU_ib(ji,jj) = 0._wp 370 zTauV_ib(ji,jj) = 0._wp 371 END DO 372 END DO 373 ENDIF 374 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 375 330 376 !------------------------------------------------------------------------------! 331 377 ! 3) Solution of the momentum equation, iterative procedure … … 391 437 ENDIF 392 438 393 ! stress at T points 394 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta) ) * z1_alph1395 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2439 ! stress at T points (zkt/=0 if landfast) 440 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 441 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 396 442 397 443 END DO … … 412 458 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) ) 413 459 414 ! stress at F points 415 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2460 ! stress at F points (zkt/=0 if landfast) 461 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 416 462 417 463 END DO … … 461 507 ! 462 508 ! !--- tau_bottom/v_ice 463 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj)) )464 zTauB = - tau_icebfr(ji,jj) / zvel509 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 510 zTauB = - zTauV_ib(ji,jj) / zvel 465 511 ! 466 512 ! !--- Coriolis at V-points (energy conserving formulation) … … 473 519 ! 474 520 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 475 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )521 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 476 522 ! 477 523 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 509 555 ! 510 556 ! !--- tau_bottom/u_ice 511 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj)) )512 zTauB = - tau_icebfr(ji,jj) / zvel557 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 558 zTauB = - zTauU_ib(ji,jj) / zvel 513 559 ! 514 560 ! !--- Coriolis at U-points (energy conserving formulation) … … 521 567 ! 522 568 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 523 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )569 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 524 570 ! 525 571 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 559 605 ! 560 606 ! !--- tau_bottom/u_ice 561 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj)) )562 zTauB = - tau_icebfr(ji,jj) / zvel607 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 608 zTauB = - zTauU_ib(ji,jj) / zvel 563 609 ! 564 610 ! !--- Coriolis at U-points (energy conserving formulation) … … 571 617 ! 572 618 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 573 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )619 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 574 620 ! 575 621 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 607 653 ! 608 654 ! !--- tau_bottom/v_ice 609 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj)) )610 z tauB = - tau_icebfr(ji,jj) / zvel655 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 656 zTauB = - zTauV_ib(ji,jj) / zvel 611 657 ! 612 658 ! !--- Coriolis at v-points (energy conserving formulation) … … 619 665 ! 620 666 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 621 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )667 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 622 668 ! 623 669 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
Note: See TracChangeset
for help on using the changeset viewer.
