Changeset 14072 for NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 (modified) (42 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r14005 r14072 6 6 !! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code 7 7 !! 3.0 ! 2008-03 (M. Vancoppenolle) adaptation to new model 8 !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 8 !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 9 9 !! 3.3 ! 2009-05 (G.Garric) addition of the evp case 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 11 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + mEVP (Bouillon 2013) … … 28 28 USE icevar ! ice_var_sshdyn 29 29 USE icedyn_rdgrft ! sea-ice: ice strength 30 USE bdy_oce , ONLY : ln_bdy 31 USE bdyice 30 USE bdy_oce , ONLY : ln_bdy 31 USE bdyice 32 32 #if defined key_agrif 33 33 USE agrif_ice_interp … … 69 69 !! 70 70 !! ** purpose : determines sea ice drift from wind stress, ice-ocean 71 !! stress and sea-surface slope. Ice-ice interaction is described by 72 !! a non-linear elasto-viscous-plastic (EVP) law including shear 73 !! strength and a bulk rheology (Hunke and Dukowicz, 2002). 71 !! stress and sea-surface slope. Ice-ice interaction is described by 72 !! a non-linear elasto-viscous-plastic (EVP) law including shear 73 !! strength and a bulk rheology (Hunke and Dukowicz, 2002). 74 74 !! 75 75 !! The points in the C-grid look like this, dear reader … … 79 79 !! | 80 80 !! (ji-1,jj) | (ji,jj) 81 !! --------- 81 !! --------- 82 82 !! | | 83 83 !! | (ji,jj) |------(ji,jj) 84 84 !! | | 85 !! --------- 85 !! --------- 86 86 !! (ji-1,jj-1) (ji,jj-1) 87 87 !! … … 90 90 !! snow total volume (vt_s) per unit area 91 91 !! 92 !! ** Action : - compute u_ice, v_ice : the components of the 92 !! ** Action : - compute u_ice, v_ice : the components of the 93 93 !! sea-ice velocity vector 94 94 !! - compute delta_i, shear_i, divu_i, which are inputs … … 96 96 !! 97 97 !! ** Steps : 0) compute mask at F point 98 !! 1) Compute ice snow mass, ice strength 98 !! 1) Compute ice snow mass, ice strength 99 99 !! 2) Compute wind, oceanic stresses, mass terms and 100 100 !! coriolis terms of the momentum equation … … 152 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 153 153 ! ! ocean surface (ssh_m) if ice is not embedded 154 ! ! ice bottom surface if ice is embedded 154 ! ! ice bottom surface if ice is embedded 155 155 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 156 156 REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points … … 172 172 !! --- diags 173 173 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 175 175 !! --- SIMIP diags 176 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 179 179 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 180 180 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 181 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 181 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 182 182 !!------------------------------------------------------------------- 183 183 … … 229 229 ! 1) define some variables and initialize arrays 230 230 !------------------------------------------------------------------------------! 231 zrhoco = rho0 * rn_cio 231 zrhoco = rho0 * rn_cio 232 232 233 233 ! ecc2: square of yield ellipse eccenticrity … … 248 248 ENDIF 249 249 z1_dtevp = 1._wp / zdtevp 250 251 ! Initialise stress tensor 252 zs1 (:,:) = pstress1_i (:,:) 250 251 ! Initialise stress tensor 252 zs1 (:,:) = pstress1_i (:,:) 253 253 zs2 (:,:) = pstress2_i (:,:) 254 254 zs12(:,:) = pstress12_i(:,:) … … 292 292 ! dt/m at T points (for alpha and beta coefficients) 293 293 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 294 294 295 295 ! m/dt 296 296 zmU_t(ji,jj) = zmassU * z1_dtevp 297 297 zmV_t(ji,jj) = zmassV * z1_dtevp 298 298 299 299 ! Drag ice-atm. 300 300 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) … … 350 350 ! ! ==================== ! 351 351 DO jter = 1 , nn_nevp ! loop over jter ! 352 ! ! ==================== ! 352 ! ! ==================== ! 353 353 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 354 354 ! … … 377 377 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 378 378 & ) * 0.25_wp * r1_e1e2t(ji,jj) 379 379 380 380 ! divergence at T points 381 381 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & … … 383 383 & ) * r1_e1e2t(ji,jj) 384 384 zdiv2 = zdiv * zdiv 385 385 386 386 ! tension at T points 387 387 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & … … 389 389 & ) * r1_e1e2t(ji,jj) 390 390 zdt2 = zdt * zdt 391 391 392 392 ! delta at T points 393 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 393 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 394 394 395 395 END_2D … … 407 407 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 408 408 & ) * r1_e1e2t(ji,jj) 409 409 410 410 ! tension at T points (duplication to avoid communications) 411 411 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 412 412 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 413 413 & ) * r1_e1e2t(ji,jj) 414 414 415 415 ! alpha for aEVP 416 416 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m … … 427 427 ! zalph2 = zalph1 428 428 ENDIF 429 429 430 430 ! stress at T points (zkt/=0 if landfast) 431 431 zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 432 432 zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 433 433 434 434 END_2D 435 435 … … 440 440 END_2D 441 441 ENDIF 442 442 443 443 DO_2D( 1, 0, 1, 0 ) 444 444 … … 451 451 ! zalph2 = zalph2 - 1._wp 452 452 ENDIF 453 453 454 454 ! P/delta at F points 455 455 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) ) 456 456 457 457 ! stress at F points (zkt/=0 if landfast) 458 458 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 … … 519 519 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 520 520 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 521 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 521 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 522 522 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 523 523 & ) / ( zbetav + 1._wp ) & … … 574 574 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 575 575 & ) / ( zbetau + 1._wp ) & 576 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 577 & ) * zmsk00x(ji,jj) 578 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 579 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 580 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 581 & ) / MAX( zepsi, zmU_t(ji,jj) + 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_lf_relax ) & ! static friction => slow decrease to v=0 583 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 584 & ) * zmsk00x(ji,jj) 585 ENDIF 586 END_2D 587 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 588 ! 589 #if defined key_agrif 590 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 591 CALL agrif_interp_ice( 'U' ) 592 #endif 593 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 594 ! 595 ELSE ! odd iterations 596 ! 597 DO_2D( 0, 0, 0, 0 ) 598 ! !--- tau_io/(u_oce - u_ice) 599 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 600 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 601 ! !--- Ocean-to-Ice stress 602 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 603 ! 604 ! !--- tau_bottom/u_ice 605 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 606 zTauB = ztaux_base(ji,jj) / zvel 607 ! !--- OceanBottom-to-Ice stress 608 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 609 ! 610 ! !--- Coriolis at U-points (energy conserving formulation) 611 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 612 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 613 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 614 ! 615 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 616 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 617 ! 618 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 619 ! 1 = sliding friction : TauB < RHS 620 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 621 ! 622 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 623 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 624 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 625 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 626 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 627 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 628 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 629 & ) / ( zbetau + 1._wp ) & 630 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 576 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 631 577 & ) * zmsk00x(ji,jj) 632 578 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 647 593 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 648 594 ! 595 ELSE ! odd iterations 596 ! 597 DO_2D( 0, 0, 0, 0 ) 598 ! !--- tau_io/(u_oce - u_ice) 599 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 600 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 601 ! !--- Ocean-to-Ice stress 602 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 603 ! 604 ! !--- tau_bottom/u_ice 605 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 606 zTauB = ztaux_base(ji,jj) / zvel 607 ! !--- OceanBottom-to-Ice stress 608 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 609 ! 610 ! !--- Coriolis at U-points (energy conserving formulation) 611 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 612 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 613 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 614 ! 615 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 616 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 617 ! 618 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 619 ! 1 = sliding friction : TauB < RHS 620 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 621 ! 622 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 623 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 624 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 625 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 626 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 627 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 628 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 629 & ) / ( zbetau + 1._wp ) & 630 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 631 & ) * zmsk00x(ji,jj) 632 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 633 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 634 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 635 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 636 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 637 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 638 & ) * zmsk00x(ji,jj) 639 ENDIF 640 END_2D 641 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 642 ! 643 #if defined key_agrif 644 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 645 CALL agrif_interp_ice( 'U' ) 646 #endif 647 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 648 ! 649 649 DO_2D( 0, 0, 0, 0 ) 650 650 ! !--- tau_io/(v_oce - v_ice) … … 679 679 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 680 680 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 681 & ) / ( zbetav + 1._wp ) & 681 & ) / ( zbetav + 1._wp ) & 682 682 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 683 683 & ) * zmsk00y(ji,jj) … … 710 710 ! 711 711 !------------------------------------------------------------------------------! 712 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 712 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 713 713 !------------------------------------------------------------------------------! 714 714 DO_2D( 1, 0, 1, 0 ) … … 720 720 721 721 END_2D 722 722 723 723 DO_2D( 0, 0, 0, 0 ) ! no vector loop 724 724 725 725 ! tension**2 at T points 726 726 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & … … 730 730 731 731 zten_i(ji,jj) = zdt 732 732 733 733 ! shear**2 at T points (doc eq. A16) 734 734 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 735 735 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 736 736 & ) * 0.25_wp * r1_e1e2t(ji,jj) 737 737 738 738 ! shear at T points 739 739 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) … … 743 743 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 744 744 & ) * r1_e1e2t(ji,jj) 745 745 746 746 ! delta at T points 747 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 747 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 748 748 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 749 749 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl … … 752 752 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 753 753 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 754 754 755 755 ! --- Store the stress tensor for the next time step --- ! 756 756 pstress1_i (:,:) = zs1 (:,:) … … 776 776 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 777 777 ENDIF 778 778 779 779 ! --- divergence, shear and strength --- ! 780 780 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence … … 786 786 ! 787 787 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 788 ! 788 ! 789 789 DO_2D( 1, 1, 1, 1 ) 790 790 791 791 ! Ice stresses 792 792 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 793 793 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 794 794 ! I know, this can be confusing... 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 796 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 797 797 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 798 798 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 799 799 800 800 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 801 801 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 802 802 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 803 804 END_2D 803 804 END_2D 805 805 ! 806 806 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 807 807 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 808 808 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 809 809 810 810 DEALLOCATE ( zsig_I, zsig_II ) 811 811 812 812 ENDIF 813 813 … … 818 818 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 819 819 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 822 822 DO_2D( 1, 1, 1, 1 ) 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 825 825 ! and **deformations** at current iterates 826 826 ! following Lemieux & Dupont (2020) … … 829 829 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 830 830 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 831 831 832 832 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 833 833 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 834 834 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 835 835 836 836 ! Normalized principal stresses (used to display the ellipse) 837 837 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 838 838 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 839 839 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 840 END_2D 841 ! 842 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 843 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 840 END_2D 841 ! 842 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 843 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 844 844 845 845 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 846 846 847 847 ENDIF 848 848 … … 889 889 890 890 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) 891 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 891 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 892 892 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 893 893 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport … … 911 911 ENDIF 912 912 ENDIF 913 ENDIF 913 ENDIF 914 914 ! 915 915 DEALLOCATE( zmsk00, zmsk15 ) … … 921 921 !!---------------------------------------------------------------------- 922 922 !! *** ROUTINE rhg_cvg *** 923 !! 923 !! 924 924 !! ** Purpose : check convergence of oce rheology 925 925 !! … … 929 929 !! This routine is called every sub-iteration, so it is cpu expensive 930 930 !! 931 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 931 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 932 932 !!---------------------------------------------------------------------- 933 933 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index … … 936 936 INTEGER :: it, idtime, istatus 937 937 INTEGER :: ji, jj ! dummy loop indices 938 REAL(wp) :: zresm ! local real 938 REAL(wp) :: zresm ! local real 939 939 CHARACTER(len=20) :: clname 940 940 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence … … 963 963 ! time 964 964 it = ( kt - 1 ) * kitermax + kiter 965 965 966 966 ! convergence 967 967 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) … … 982 982 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 983 983 ENDIF 984 984 985 985 END SUBROUTINE rhg_cvg 986 986 … … 989 989 !!--------------------------------------------------------------------- 990 990 !! *** ROUTINE rhg_evp_rst *** 991 !! 991 !! 992 992 !! ** Purpose : Read or write RHG file in restart file 993 993 !! … … 1041 1041 END SUBROUTINE rhg_evp_rst 1042 1042 1043 1043 1044 1044 #else 1045 1045 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
