Changeset 13746 for NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90
- Timestamp:
- 2020-11-08T22:25:28+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90
r13574 r13746 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 49 50 !! * Substitutions 50 51 # include "vectopt_loop_substitute.h90" 52 53 !! for convergence tests 54 INTEGER :: ncvgid ! netcdf file id 55 INTEGER :: nvarid ! netcdf variable id 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 51 57 !!---------------------------------------------------------------------- 52 58 !! NEMO/ICE 4.0 , NEMO Consortium (2018) 53 !! $Id: icedyn_rhg_evp.F90 1 1536 2019-09-11 13:54:18Z smasson$59 !! $Id: icedyn_rhg_evp.F90 13662 2020-10-22 18:49:56Z clem $ 54 60 !! Software governed by the CeCILL license (see ./LICENSE) 55 61 !!---------------------------------------------------------------------- … … 119 125 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 120 126 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 127 REAl(wp) :: zbetau, zbetav 121 128 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 122 REAL(wp) :: z delta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars129 REAL(wp) :: zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 123 130 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 124 131 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 125 132 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 126 133 ! 127 REAL(wp) :: zresm ! Maximal error on ice velocity128 134 REAL(wp) :: zintb, zintn ! dummy argument 129 135 REAL(wp) :: zfac_x, zfac_y 130 136 REAL(wp) :: zshear, zdum1, zdum2 131 REAL(wp) :: invw ! for test case 132 ! 133 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 137 REAL(wp) :: zinvw ! for test case 138 139 ! 140 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 141 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 134 142 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 135 143 ! … … 138 146 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 139 147 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 140 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 148 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 141 149 ! 142 150 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 143 151 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 144 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence145 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 146 153 ! ! ocean surface (ssh_m) if ice is not embedded … … 156 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 157 164 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 158 REAL(wp), DIMENSION(jpi,jpj) :: zfmask , zwf! mask at F points for the ice165 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 159 166 160 167 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 161 168 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 162 169 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 170 !! --- check convergence 171 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 163 172 !! --- diags 164 REAL(wp) , DIMENSION(jpi,jpj) :: zmsk00165 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig 1, zsig2, zsig3173 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 166 175 !! --- SIMIP diags 167 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 175 184 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 176 185 ! 177 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 186 ! for diagnostics and convergence tests 187 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 188 DO jj = 1, jpj 189 DO ji = 1, jpi 190 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 191 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 192 END DO 193 END DO 194 ! 195 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 178 196 !------------------------------------------------------------------------------! 179 197 ! 0) mask at F points for the ice … … 188 206 189 207 ! Lateral boundary conditions on velocity (modify zfmask) 190 zwf(:,:) = zfmask(:,:)191 208 DO jj = 2, jpjm1 192 209 DO ji = fs_2, fs_jpim1 ! vector opt. 193 210 IF( zfmask(ji,jj) == 0._wp ) THEN 194 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 211 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 212 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 195 213 ENDIF 196 214 END DO … … 198 216 DO jj = 2, jpjm1 199 217 IF( zfmask(1,jj) == 0._wp ) THEN 200 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )218 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 201 219 ENDIF 202 220 IF( zfmask(jpi,jj) == 0._wp ) THEN 203 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )204 221 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 222 ENDIF 205 223 END DO 206 224 DO ji = 2, jpim1 207 225 IF( zfmask(ji,1) == 0._wp ) THEN 208 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )226 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 209 227 ENDIF 210 228 IF( zfmask(ji,jpj) == 0._wp ) THEN 211 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )229 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 212 230 ENDIF 213 231 END DO … … 219 237 zrhoco = rau0 * rn_cio 220 238 !extra code for test case boundary conditions 221 invw=1._wp/(zrhoco*0.5_wp)239 zinvw=1._wp/(zrhoco*0.5_wp) 222 240 223 241 ! ecc2: square of yield ellipse eccenticrity … … 225 243 z1_ecc2 = 1._wp / ecc2 226 244 227 ! Time step for subcycling228 zdtevp = rdt_ice / REAL( nn_nevp )229 z1_dtevp = 1._wp / zdtevp230 231 245 ! alpha parameters (Bouillon 2009) 232 246 IF( .NOT. ln_aEVP ) THEN 233 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 247 zdtevp = rdt_ice / REAL( nn_nevp ) 248 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 234 249 zalph2 = zalph1 * z1_ecc2 235 250 236 251 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 237 252 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 253 ELSE 254 zdtevp = rdt_ice 255 ! zalpha parameters set later on adaptatively 238 256 ENDIF 257 z1_dtevp = 1._wp / zdtevp 239 258 240 259 ! Initialise stress tensor … … 247 266 248 267 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 249 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile268 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 250 269 ELSE ; zkt = 0._wp 251 270 ENDIF … … 318 337 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) 319 338 ! ice-bottom stress at U points 320 zvCr = zaU(ji,jj) * rn_ depfra * hu_n(ji,jj)321 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )339 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 340 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 322 341 ! ice-bottom stress at V points 323 zvCr = zaV(ji,jj) * rn_ depfra * hv_n(ji,jj)324 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )342 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 343 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 325 344 ! ice_bottom stress at T points 326 zvCr = at_i(ji,jj) * rn_ depfra * ht_n(ji,jj)327 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )345 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 346 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 328 347 END DO 329 348 END DO … … 348 367 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 349 368 ! 350 !!$ IF(ln_ctl) THEN ! Convergence test 351 !!$ DO jj = 1, jpjm1 352 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 353 !!$ zv_ice(:,jj) = v_ice(:,jj) 354 !!$ END DO 355 !!$ ENDIF 369 ! convergence test 370 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 371 DO jj = 1, jpj 372 DO ji = 1, jpi 373 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 374 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 375 END DO 376 END DO 377 ENDIF 356 378 357 379 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 358 DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points380 DO jj = 1, jpjm1 359 381 DO ji = 1, jpim1 360 382 … … 366 388 END DO 367 389 END DO 368 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. ) 369 370 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 371 DO ji = 2, jpi ! no vector loop 390 391 DO jj = 2, jpjm1 392 DO ji = 2, jpim1 ! no vector loop 372 393 373 394 ! shear**2 at T points (doc eq. A16) … … 389 410 390 411 ! delta at T points 391 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 392 393 ! P/delta at T points 394 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 395 396 ! alpha & beta for aEVP 412 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 413 414 END DO 415 END DO 416 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1._wp ) 417 418 ! P/delta at T points 419 DO jj = 1, jpj 420 DO ji = 1, jpi 421 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 422 END DO 423 END DO 424 425 DO jj = 2, jpj ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 426 DO ji = 2, jpi ! no vector loop 427 428 ! divergence at T points (duplication to avoid communications) 429 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 430 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 431 & ) * r1_e1e2t(ji,jj) 432 433 ! tension at T points (duplication to avoid communications) 434 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) & 435 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 436 & ) * r1_e1e2t(ji,jj) 437 438 ! alpha for aEVP 397 439 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 398 440 ! alpha = beta = sqrt(4*gamma) … … 402 444 zalph2 = zalph1 403 445 z1_alph2 = z1_alph1 446 ! explicit: 447 ! z1_alph1 = 1._wp / zalph1 448 ! z1_alph2 = 1._wp / zalph1 449 ! zalph1 = zalph1 - 1._wp 450 ! zalph2 = zalph1 404 451 ENDIF 405 452 406 453 ! stress at T points (zkt/=0 if landfast) 407 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta *(1._wp - zkt) ) ) * z1_alph1408 zs2(ji,jj) = ( zs2(ji,jj) *zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2454 zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 455 zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 409 456 410 457 END DO 411 458 END DO 412 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 413 459 460 ! Save beta at T-points for further computations 461 IF( ln_aEVP ) THEN 462 DO jj = 1, jpj 463 DO ji = 1, jpi 464 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 465 END DO 466 END DO 467 ENDIF 468 414 469 DO jj = 1, jpjm1 415 470 DO ji = 1, jpim1 416 471 417 ! alpha & betafor aEVP472 ! alpha for aEVP 418 473 IF( ln_aEVP ) THEN 419 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj)) )474 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 420 475 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 421 zbeta(ji,jj) = zalph2 476 ! explicit: 477 ! z1_alph2 = 1._wp / zalph2 478 ! zalph2 = zalph2 - 1._wp 422 479 ENDIF 423 480 … … 489 546 ! 490 547 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 491 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 492 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 493 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 494 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 495 & ) * 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 548 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 549 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 550 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 551 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 552 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 553 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 554 & ) / ( zbetav + 1._wp ) & 555 & ) * 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 496 556 & ) * zmsk00y(ji,jj) 497 557 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 498 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity499 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)500 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast501 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0502 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin503 & ) 558 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 559 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 560 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 561 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 562 & ) * 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 563 & ) * zmsk00y(ji,jj) 504 564 ENDIF 505 565 !extra code for test case boundary conditions 506 566 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 507 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))567 v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 508 568 END IF 509 569 END DO … … 544 604 ! 545 605 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 546 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 547 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 548 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 549 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 550 & ) * 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 606 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 607 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 608 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 609 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 610 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 611 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 612 & ) / ( zbetau + 1._wp ) & 613 & ) * 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 551 614 & ) * zmsk00x(ji,jj) 552 615 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 553 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity554 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)555 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast556 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0557 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin558 & 616 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 617 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 618 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 619 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 620 & ) * 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 621 & ) * zmsk00x(ji,jj) 559 622 ENDIF 560 623 !extra code for test case boundary conditions 561 624 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 562 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))625 u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 563 626 END IF 564 627 END DO … … 601 664 ! 602 665 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 603 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 604 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 605 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 606 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 607 & ) * 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 666 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 667 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 668 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 669 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 670 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 671 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 672 & ) / ( zbetau + 1._wp ) & 673 & ) * 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 608 674 & ) * zmsk00x(ji,jj) 609 675 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 610 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity611 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)612 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast613 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0614 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin615 & 676 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 677 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 678 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 679 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 680 & ) * 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 681 & ) * zmsk00x(ji,jj) 616 682 ENDIF 617 683 !extra code for test case boundary conditions 618 684 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 619 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))685 u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 620 686 END IF 621 687 END DO … … 656 722 ! 657 723 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 658 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 659 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 660 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 661 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 662 & ) * 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 724 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 725 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 726 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 727 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 728 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 729 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 730 & ) / ( zbetav + 1._wp ) & 731 & ) * 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 663 732 & ) * zmsk00y(ji,jj) 664 733 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 665 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity666 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)667 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast668 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0669 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin670 & 734 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 735 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 736 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 737 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 738 & ) * 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 739 & ) * zmsk00y(ji,jj) 671 740 ENDIF 672 741 !extra code for test case boundary conditions 673 742 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 674 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))743 v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 675 744 END IF 676 745 END DO … … 686 755 ENDIF 687 756 688 !!$ IF(ln_ctl) THEN ! Convergence test 689 !!$ DO jj = 2 , jpjm1 690 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 691 !!$ END DO 692 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 693 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 694 !!$ ENDIF 757 ! convergence test 758 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 695 759 ! 696 760 ! ! ==================== ! 697 761 END DO ! end loop over jter ! 698 762 ! ! ==================== ! 763 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 699 764 ! 700 765 !------------------------------------------------------------------------------! … … 721 786 zdt2 = zdt * zdt 722 787 788 zten_i(ji,jj) = zdt 789 723 790 ! shear**2 at T points (doc eq. A16) 724 791 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 735 802 736 803 ! delta at T points 737 z delta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )738 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0739 pdelta_i(ji,jj) = z delta + rn_creepl * rswitch804 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 805 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 806 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 740 807 741 808 END DO 742 809 END DO 743 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 810 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 811 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. ) 744 812 745 813 ! --- Store the stress tensor for the next time step --- ! 746 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )747 814 pstress1_i (:,:) = zs1 (:,:) 748 815 pstress2_i (:,:) = zs2 (:,:) … … 753 820 ! 5) diagnostics 754 821 !------------------------------------------------------------------------------! 755 DO jj = 1, jpj756 DO ji = 1, jpi757 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice758 END DO759 END DO760 761 822 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 762 823 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 777 838 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 778 839 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 840 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 779 841 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 780 842 781 ! --- stress tensor--- !782 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN783 ! 784 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )843 ! --- Stress tensor invariants (SIMIP diags) --- ! 844 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 845 ! 846 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 785 847 ! 786 DO jj = 2, jpjm1 787 DO ji = 2, jpim1 788 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 789 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 790 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 791 792 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 793 794 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 795 796 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 797 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 798 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 799 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 800 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 801 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 802 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 803 END DO 804 END DO 805 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 806 ! 807 CALL iom_put( 'isig1' , zsig1 ) 808 CALL iom_put( 'isig2' , zsig2 ) 809 CALL iom_put( 'isig3' , zsig3 ) 810 ! 811 ! Stress tensor invariants (normal and shear stress N/m) 812 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 813 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 814 815 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 848 DO jj = 1, jpj 849 DO ji = 1, jpi 850 851 ! Ice stresses 852 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 853 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 854 ! I know, this can be confusing... 855 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 856 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 857 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 858 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 859 860 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 861 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 862 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 863 864 END DO 865 END DO 866 ! 867 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 868 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 869 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 870 871 DEALLOCATE ( zsig_I, zsig_II ) 872 816 873 ENDIF 874 875 ! --- Normalized stress tensor principal components --- ! 876 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 877 ! Recommendation 1 : we use ice strength, not replacement pressure 878 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 879 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 880 ! 881 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 882 ! 883 DO jj = 1, jpj 884 DO ji = 1, jpi 885 886 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 887 ! and **deformations** at current iterates 888 ! following Lemieux & Dupont (2020) 889 zfac = zp_delt(ji,jj) 890 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 891 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 892 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 893 894 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 895 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 896 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 817 897 898 ! Normalized principal stresses (used to display the ellipse) 899 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 900 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 901 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 902 END DO 903 END DO 904 ! 905 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 906 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 907 908 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 909 910 ENDIF 911 818 912 ! --- SIMIP --- ! 819 913 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 871 965 ENDIF 872 966 ! 967 ! --- convergence tests --- ! 968 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 969 IF( iom_use('uice_cvg') ) THEN 970 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 971 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 972 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 973 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 974 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 975 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 976 ENDIF 977 ENDIF 978 ENDIF 979 ! 980 DEALLOCATE( zmsk00, zmsk15 ) 981 ! 873 982 END SUBROUTINE ice_dyn_rhg_evp 983 984 985 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 986 !!---------------------------------------------------------------------- 987 !! *** ROUTINE rhg_cvg *** 988 !! 989 !! ** Purpose : check convergence of oce rheology 990 !! 991 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 992 !! during the sub timestepping of rheology so as: 993 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 994 !! This routine is called every sub-iteration, so it is cpu expensive 995 !! 996 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 997 !!---------------------------------------------------------------------- 998 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 999 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 1000 !! 1001 INTEGER :: it, idtime, istatus 1002 INTEGER :: ji, jj ! dummy loop indices 1003 REAL(wp) :: zresm ! local real 1004 CHARACTER(len=20) :: clname 1005 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 1006 !!---------------------------------------------------------------------- 1007 1008 ! create file 1009 IF( kt == nit000 .AND. kiter == 1 ) THEN 1010 ! 1011 IF( lwp ) THEN 1012 WRITE(numout,*) 1013 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 1014 WRITE(numout,*) '~~~~~~~' 1015 ENDIF 1016 ! 1017 IF( lwm ) THEN 1018 clname = 'ice_cvg.nc' 1019 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 1020 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 1021 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 1022 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 1023 istatus = NF90_ENDDEF(ncvgid) 1024 ENDIF 1025 ! 1026 ENDIF 1027 1028 ! time 1029 it = ( kt - 1 ) * kitermax + kiter 1030 1031 ! convergence 1032 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 1033 zresm = 0._wp 1034 ELSE 1035 DO jj = 1, jpj 1036 DO ji = 1, jpi 1037 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1038 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 1039 END DO 1040 END DO 1041 1042 ! cut of the boundary of the box (forced velocities) 1043 IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN 1044 zres(ji,jj) = 0._wp 1045 END IF 1046 1047 zresm = MAXVAL( zres ) 1048 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1049 ENDIF 1050 1051 IF( lwm ) THEN 1052 ! write variables 1053 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1054 ! close file 1055 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 1056 ENDIF 1057 1058 END SUBROUTINE rhg_cvg 874 1059 875 1060 … … 929 1114 END SUBROUTINE rhg_evp_rst 930 1115 1116 931 1117 #else 932 1118 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.