Changeset 10419 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2018-12-19T20:46:30+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90
r10402 r10419 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 … … 144 147 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 145 148 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence 146 REAL(wp), DIMENSION(jpi,jpj) :: zssh _lead_m! array used for the calculation of ice surface slope:149 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 147 150 ! ! ocean surface (ssh_m) if ice is not embedded 148 ! ! ice topsurface if ice is embedded151 ! ! ice bottom surface if ice is embedded 149 152 REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array 150 153 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array … … 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 !------------------------------------------------------------------------------! 261 268 ! 2) Wind / ocean stress, mass terms, coriolis terms 262 269 !------------------------------------------------------------------------------! 263 264 ! == embedded sea ice: compute representative ice top surface ==!265 ! == non-embedded sea ice: use ocean surface for slope calculation ==!266 zssh _lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)270 ! sea surface height 271 ! embedded sea ice: compute representative ice top surface 272 ! non-embedded sea ice: use ocean surface for slope calculation 273 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 267 274 268 275 DO jj = 2, jpjm1 … … 302 309 303 310 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 304 zspgU(ji,jj) = - zmassU * grav * ( zssh _lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj)305 zspgV(ji,jj) = - zmassV * grav * ( zssh _lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj)311 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 312 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 306 313 307 314 ! masks … … 317 324 CALL lbc_lnk_multi( 'icedyn_rhg_evp', 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( 'icedyn_rhg_evp', 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 … … 382 428 ENDIF 383 429 384 ! stress at T points 385 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta) ) * z1_alph1386 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2430 ! stress at T points (zkt/=0 if landfast) 431 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 432 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 387 433 388 434 END DO … … 403 449 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) ) 404 450 405 ! stress at F points 406 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2451 ! stress at F points (zkt/=0 if landfast) 452 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 407 453 408 454 END DO … … 452 498 ! 453 499 ! !--- tau_bottom/v_ice 454 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj)) )455 zTauB = - tau_icebfr(ji,jj) / zvel500 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 501 zTauB = - zTauV_ib(ji,jj) / zvel 456 502 ! 457 503 ! !--- Coriolis at V-points (energy conserving formulation) … … 464 510 ! 465 511 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 466 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )512 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 467 513 ! 468 514 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 500 546 ! 501 547 ! !--- tau_bottom/u_ice 502 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj)) )503 zTauB = - tau_icebfr(ji,jj) / zvel548 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 549 zTauB = - zTauU_ib(ji,jj) / zvel 504 550 ! 505 551 ! !--- Coriolis at U-points (energy conserving formulation) … … 512 558 ! 513 559 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 514 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )560 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 515 561 ! 516 562 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 550 596 ! 551 597 ! !--- tau_bottom/u_ice 552 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj)) )553 zTauB = - tau_icebfr(ji,jj) / zvel598 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 599 zTauB = - zTauU_ib(ji,jj) / zvel 554 600 ! 555 601 ! !--- Coriolis at U-points (energy conserving formulation) … … 562 608 ! 563 609 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 564 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )610 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 565 611 ! 566 612 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) … … 598 644 ! 599 645 ! !--- tau_bottom/v_ice 600 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj)) )601 z tauB = - tau_icebfr(ji,jj) / zvel646 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 647 zTauB = - zTauV_ib(ji,jj) / zvel 602 648 ! 603 649 ! !--- Coriolis at v-points (energy conserving formulation) … … 610 656 ! 611 657 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 612 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )658 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 613 659 ! 614 660 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
Note: See TracChangeset
for help on using the changeset viewer.
