Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icedyn_rhg_evp.F90
r13612 r14789 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 … … 199 199 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 200 200 END_2D 201 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp 201 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp) 202 202 203 203 ! Lateral boundary conditions on velocity (modify zfmask) … … 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) … … 316 316 317 317 END_2D 318 CALL lbc_lnk _multi( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )318 CALL lbc_lnk( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp ) 319 319 ! 320 320 ! !== Landfast ice parameterization ==! … … 326 326 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) 327 327 ! ice-bottom stress at U points 328 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 328 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 329 329 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 330 330 ! ice-bottom stress at V points 331 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 331 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 332 332 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 333 333 ! ice_bottom stress at T points 334 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 334 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 335 335 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 336 336 END_2D … … 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 750 750 751 751 END_2D 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 & 754 752 CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 753 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 754 755 755 ! --- Store the stress tensor for the next time step --- ! 756 756 pstress1_i (:,:) = zs1 (:,:) … … 766 766 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 767 767 ! 768 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, & 769 & ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 768 CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, & 769 & ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, & 770 & ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 770 771 ! 771 772 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 776 777 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 777 778 ENDIF 778 779 779 780 ! --- divergence, shear and strength --- ! 780 781 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence … … 786 787 ! 787 788 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 788 ! 789 ! 789 790 DO_2D( 1, 1, 1, 1 ) 790 791 791 792 ! Ice stresses 792 793 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 793 794 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 794 795 ! I know, this can be confusing... 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 797 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 797 798 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 798 799 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 799 800 800 801 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 801 802 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 802 803 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 803 804 END_2D 804 805 END_2D 805 806 ! 806 807 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 807 808 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 808 809 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 809 810 810 811 DEALLOCATE ( zsig_I, zsig_II ) 811 812 812 813 ENDIF 813 814 … … 818 819 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 819 820 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 821 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 822 ! 822 823 DO_2D( 1, 1, 1, 1 ) 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 824 825 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 825 826 ! and **deformations** at current iterates 826 827 ! following Lemieux & Dupont (2020) … … 829 830 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 830 831 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 831 832 832 833 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 833 834 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 834 835 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 835 836 836 837 ! Normalized principal stresses (used to display the ellipse) 837 838 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 838 839 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 839 840 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 ) 841 END_2D 842 ! 843 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 844 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 844 845 845 846 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 846 847 847 848 ENDIF 848 849 … … 851 852 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 852 853 ! 853 CALL lbc_lnk _multi( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &854 & 854 CALL lbc_lnk( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, & 855 & zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 855 856 856 857 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) … … 884 885 END_2D 885 886 886 CALL lbc_lnk _multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &887 & 888 & 887 CALL lbc_lnk( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, & 888 & zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, & 889 & zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp ) 889 890 890 891 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 892 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 892 893 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 893 894 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport … … 911 912 ENDIF 912 913 ENDIF 913 ENDIF 914 ENDIF 914 915 ! 915 916 DEALLOCATE( zmsk00, zmsk15 ) … … 921 922 !!---------------------------------------------------------------------- 922 923 !! *** ROUTINE rhg_cvg *** 923 !! 924 !! 924 925 !! ** Purpose : check convergence of oce rheology 925 926 !! … … 929 930 !! This routine is called every sub-iteration, so it is cpu expensive 930 931 !! 931 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 932 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 932 933 !!---------------------------------------------------------------------- 933 934 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index … … 936 937 INTEGER :: it, idtime, istatus 937 938 INTEGER :: ji, jj ! dummy loop indices 938 REAL(wp) :: zresm ! local real 939 REAL(wp) :: zresm ! local real 939 940 CHARACTER(len=20) :: clname 940 941 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence … … 963 964 ! time 964 965 it = ( kt - 1 ) * kitermax + kiter 965 966 966 967 ! convergence 967 968 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) … … 982 983 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 983 984 ENDIF 984 985 985 986 END SUBROUTINE rhg_cvg 986 987 … … 989 990 !!--------------------------------------------------------------------- 990 991 !! *** ROUTINE rhg_evp_rst *** 991 !! 992 !! 992 993 !! ** Purpose : Read or write RHG file in restart file 993 994 !! … … 1033 1034 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 1034 1035 ! 1035 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i 1036 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i 1036 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 1037 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 1037 1038 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 1038 1039 ! … … 1041 1042 END SUBROUTINE rhg_evp_rst 1042 1043 1043 1044 1044 1045 #else 1045 1046 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
