Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icedyn_rhg_evp.F90
r13295 r14037 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 50 51 # include "do_loop_substitute.h90" 51 52 # include "domzgr_substitute.h90" 53 54 !! for convergence tests 55 INTEGER :: ncvgid ! netcdf file id 56 INTEGER :: nvarid ! netcdf variable id 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 52 58 !!---------------------------------------------------------------------- 53 59 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 121 127 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 122 128 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 129 REAl(wp) :: zbetau, zbetav 123 130 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 124 REAL(wp) :: z delta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars131 REAL(wp) :: zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 125 132 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 126 133 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 127 134 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 128 135 ! 129 REAL(wp) :: zresm ! Maximal error on ice velocity130 136 REAL(wp) :: zintb, zintn ! dummy argument 131 137 REAL(wp) :: zfac_x, zfac_y 132 138 REAL(wp) :: zshear, zdum1, zdum2 133 139 ! 134 REAL(wp), DIMENSION(jpi,jpj) :: z p_delt !P/delta at T points140 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 135 141 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 136 142 ! … … 139 145 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 140 146 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 141 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 147 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 142 148 ! 143 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 150 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 144 151 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 145 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence146 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 147 153 ! ! ocean surface (ssh_m) if ice is not embedded … … 157 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 158 164 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 159 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 160 166 161 167 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 162 168 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 163 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 164 172 !! --- diags 165 REAL(wp) , DIMENSION(jpi,jpj) :: zmsk00166 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 167 175 !! --- SIMIP diags 168 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 176 184 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 177 185 ! 178 !!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_2D( 1, 1, 1, 1 ) 189 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 190 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 191 END_2D 192 ! 193 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 179 194 !------------------------------------------------------------------------------! 180 195 ! 0) mask at F points for the ice … … 184 199 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 185 200 END_2D 186 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp 201 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp) 187 202 188 203 ! Lateral boundary conditions on velocity (modify zfmask) 189 zwf(:,:) = zfmask(:,:)190 204 DO_2D( 0, 0, 0, 0 ) 191 205 IF( zfmask(ji,jj) == 0._wp ) THEN 192 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) ) ) 206 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 207 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 193 208 ENDIF 194 209 END_2D 195 210 DO jj = 2, jpjm1 196 211 IF( zfmask(1,jj) == 0._wp ) THEN 197 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )212 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 198 213 ENDIF 199 214 IF( zfmask(jpi,jj) == 0._wp ) THEN 200 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )201 215 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 216 ENDIF 202 217 END DO 203 218 DO ji = 2, jpim1 204 219 IF( zfmask(ji,1) == 0._wp ) THEN 205 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )220 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 206 221 ENDIF 207 222 IF( zfmask(ji,jpj) == 0._wp ) THEN 208 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )223 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 209 224 ENDIF 210 225 END DO … … 220 235 z1_ecc2 = 1._wp / ecc2 221 236 222 ! Time step for subcycling223 zdtevp = rDt_ice / REAL( nn_nevp )224 z1_dtevp = 1._wp / zdtevp225 226 237 ! alpha parameters (Bouillon 2009) 227 238 IF( .NOT. ln_aEVP ) THEN 228 zalph1 = ( 2._wp * rn_relast * rDt_ice ) * z1_dtevp 239 zdtevp = rDt_ice / REAL( nn_nevp ) 240 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 229 241 zalph2 = zalph1 * z1_ecc2 230 242 231 243 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 232 244 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 233 ENDIF 245 ELSE 246 zdtevp = rdt_ice 247 ! zalpha parameters set later on adaptatively 248 ENDIF 249 z1_dtevp = 1._wp / zdtevp 234 250 235 251 ! Initialise stress tensor … … 242 258 243 259 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 244 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile260 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 245 261 ELSE ; zkt = 0._wp 246 262 ENDIF … … 310 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) 311 327 ! ice-bottom stress at U points 312 zvCr = zaU(ji,jj) * rn_ depfra * hu(ji,jj,Kmm)313 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )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 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 314 330 ! ice-bottom stress at V points 315 zvCr = zaV(ji,jj) * rn_ depfra * hv(ji,jj,Kmm)316 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )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 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 317 333 ! ice_bottom stress at T points 318 zvCr = at_i(ji,jj) * rn_ depfra * ht(ji,jj)319 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(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 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 320 336 END_2D 321 337 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp ) … … 337 353 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 338 354 ! 339 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 340 !!$ DO jj = 1, jpjm1 341 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 342 !!$ zv_ice(:,jj) = v_ice(:,jj) 343 !!$ END DO 344 !!$ ENDIF 355 ! convergence test 356 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 357 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 358 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 359 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 360 END_2D 361 ENDIF 345 362 346 363 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 353 370 354 371 END_2D 355 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 356 357 DO_2D( 0, 1, 0, 1 ) 372 373 DO_2D( 0, 0, 0, 0 ) 358 374 359 375 ! shear**2 at T points (doc eq. A16) … … 375 391 376 392 ! delta at T points 377 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 378 379 ! P/delta at T points 380 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 381 382 ! alpha & beta for aEVP 393 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 394 395 END_2D 396 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 397 398 ! P/delta at T points 399 DO_2D( 1, 1, 1, 1 ) 400 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 401 END_2D 402 403 DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 404 405 ! divergence at T points (duplication to avoid communications) 406 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 407 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 408 & ) * r1_e1e2t(ji,jj) 409 410 ! tension at T points (duplication to avoid communications) 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 & - ( 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 & ) * r1_e1e2t(ji,jj) 414 415 ! alpha for aEVP 383 416 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 384 417 ! alpha = beta = sqrt(4*gamma) … … 388 421 zalph2 = zalph1 389 422 z1_alph2 = z1_alph1 423 ! explicit: 424 ! z1_alph1 = 1._wp / zalph1 425 ! z1_alph2 = 1._wp / zalph1 426 ! zalph1 = zalph1 - 1._wp 427 ! zalph2 = zalph1 390 428 ENDIF 391 429 392 430 ! stress at T points (zkt/=0 if landfast) 393 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta *(1._wp - zkt) ) ) * z1_alph1394 zs2(ji,jj) = ( zs2(ji,jj) *zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2431 zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(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 395 433 396 434 END_2D 397 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 398 435 436 ! Save beta at T-points for further computations 437 IF( ln_aEVP ) THEN 438 DO_2D( 1, 1, 1, 1 ) 439 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 440 END_2D 441 ENDIF 442 399 443 DO_2D( 1, 0, 1, 0 ) 400 444 401 ! alpha & betafor aEVP445 ! alpha for aEVP 402 446 IF( ln_aEVP ) THEN 403 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj)) )447 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 404 448 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 405 zbeta(ji,jj) = zalph2 449 ! explicit: 450 ! z1_alph2 = 1._wp / zalph2 451 ! zalph2 = zalph2 - 1._wp 406 452 ENDIF 407 453 … … 469 515 ! 470 516 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 471 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 472 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 473 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 474 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 475 & ) * 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 517 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 518 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 519 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 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) & 522 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 523 & ) / ( zbetav + 1._wp ) & 524 & ) * 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 476 525 & ) * zmsk00y(ji,jj) 477 526 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 478 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity479 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)480 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast481 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0482 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin483 & ) 527 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 528 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 529 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 530 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 531 & ) * 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 532 & ) * zmsk00y(ji,jj) 484 533 ENDIF 485 534 END_2D … … 518 567 ! 519 568 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 520 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 521 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 522 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 523 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 524 & ) * 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 569 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 570 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 571 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 572 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 573 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 574 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 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 525 577 & ) * zmsk00x(ji,jj) 526 578 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 527 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity528 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)529 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast530 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0531 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin532 & 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) 533 585 ENDIF 534 586 END_2D … … 569 621 ! 570 622 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 571 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 572 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 573 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 574 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 575 & ) * 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 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 631 & ) * zmsk00x(ji,jj) 577 632 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 578 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity579 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)580 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast581 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0582 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin583 & 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) 584 639 ENDIF 585 640 END_2D … … 618 673 ! 619 674 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 620 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 621 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 622 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 623 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 624 & ) * 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 675 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 676 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 677 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 678 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 679 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 680 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 681 & ) / ( zbetav + 1._wp ) & 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 625 683 & ) * zmsk00y(ji,jj) 626 684 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 627 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity628 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)629 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast630 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0631 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin632 & 685 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 686 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 687 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 688 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 689 & ) * 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 690 & ) * zmsk00y(ji,jj) 633 691 ENDIF 634 692 END_2D … … 643 701 ENDIF 644 702 645 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 646 !!$ DO jj = 2 , jpjm1 647 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 648 !!$ END DO 649 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 650 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 651 !!$ ENDIF 703 ! convergence test 704 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 652 705 ! 653 706 ! ! ==================== ! 654 707 END DO ! end loop over jter ! 655 708 ! ! ==================== ! 709 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 656 710 ! 657 711 !------------------------------------------------------------------------------! … … 667 721 END_2D 668 722 669 DO_2D( 0, 0, 0, 0 ) 723 DO_2D( 0, 0, 0, 0 ) ! no vector loop 670 724 671 725 ! tension**2 at T points … … 674 728 & ) * r1_e1e2t(ji,jj) 675 729 zdt2 = zdt * zdt 730 731 zten_i(ji,jj) = zdt 676 732 677 733 ! shear**2 at T points (doc eq. A16) … … 689 745 690 746 ! delta at T points 691 z delta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )692 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0693 pdelta_i(ji,jj) = z delta + rn_creepl * rswitch747 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 748 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 749 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 694 750 695 751 END_2D 696 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 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 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 697 754 698 755 ! --- Store the stress tensor for the next time step --- ! 699 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )700 756 pstress1_i (:,:) = zs1 (:,:) 701 757 pstress2_i (:,:) = zs2 (:,:) … … 706 762 ! 5) diagnostics 707 763 !------------------------------------------------------------------------------! 708 DO_2D( 1, 1, 1, 1 )709 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice710 END_2D711 712 764 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 713 765 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 730 782 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 731 783 732 ! --- stress tensor--- !733 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN734 ! 735 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )784 ! --- Stress tensor invariants (SIMIP diags) --- ! 785 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 786 ! 787 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 736 788 ! 737 DO_2D( 0, 0, 0, 0 ) 738 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 739 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 740 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 741 742 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 743 744 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 745 746 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 747 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 748 !! 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 749 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 750 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 751 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 752 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 753 END_2D 754 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 755 ! 756 CALL iom_put( 'isig1' , zsig1 ) 757 CALL iom_put( 'isig2' , zsig2 ) 758 CALL iom_put( 'isig3' , zsig3 ) 759 ! 760 ! Stress tensor invariants (normal and shear stress N/m) 761 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 762 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 763 764 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 765 ENDIF 766 789 DO_2D( 1, 1, 1, 1 ) 790 791 ! Ice stresses 792 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 793 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 794 ! I know, this can be confusing... 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 797 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 798 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 799 800 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 801 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 802 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 803 804 END_2D 805 ! 806 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 807 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 808 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 809 810 DEALLOCATE ( zsig_I, zsig_II ) 811 812 ENDIF 813 814 ! --- Normalized stress tensor principal components --- ! 815 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 816 ! Recommendation 1 : we use ice strength, not replacement pressure 817 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 818 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 819 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 822 DO_2D( 1, 1, 1, 1 ) 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 825 ! and **deformations** at current iterates 826 ! following Lemieux & Dupont (2020) 827 zfac = zp_delt(ji,jj) 828 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 829 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 830 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 831 832 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 833 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 834 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 835 836 ! Normalized principal stresses (used to display the ellipse) 837 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 838 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 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 ) 844 845 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 846 847 ENDIF 848 767 849 ! --- SIMIP --- ! 768 850 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 818 900 ENDIF 819 901 ! 902 ! --- convergence tests --- ! 903 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 904 IF( iom_use('uice_cvg') ) THEN 905 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 906 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 907 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 908 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 909 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 910 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 911 ENDIF 912 ENDIF 913 ENDIF 914 ! 915 DEALLOCATE( zmsk00, zmsk15 ) 916 ! 820 917 END SUBROUTINE ice_dyn_rhg_evp 918 919 920 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 921 !!---------------------------------------------------------------------- 922 !! *** ROUTINE rhg_cvg *** 923 !! 924 !! ** Purpose : check convergence of oce rheology 925 !! 926 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 927 !! during the sub timestepping of rheology so as: 928 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 929 !! This routine is called every sub-iteration, so it is cpu expensive 930 !! 931 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 932 !!---------------------------------------------------------------------- 933 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 934 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 935 !! 936 INTEGER :: it, idtime, istatus 937 INTEGER :: ji, jj ! dummy loop indices 938 REAL(wp) :: zresm ! local real 939 CHARACTER(len=20) :: clname 940 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 941 !!---------------------------------------------------------------------- 942 943 ! create file 944 IF( kt == nit000 .AND. kiter == 1 ) THEN 945 ! 946 IF( lwp ) THEN 947 WRITE(numout,*) 948 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 949 WRITE(numout,*) '~~~~~~~' 950 ENDIF 951 ! 952 IF( lwm ) THEN 953 clname = 'ice_cvg.nc' 954 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 955 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 956 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 957 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 958 istatus = NF90_ENDDEF(ncvgid) 959 ENDIF 960 ! 961 ENDIF 962 963 ! time 964 it = ( kt - 1 ) * kitermax + kiter 965 966 ! convergence 967 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 968 zresm = 0._wp 969 ELSE 970 DO_2D( 1, 1, 1, 1 ) 971 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 972 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 973 END_2D 974 zresm = MAXVAL( zres ) 975 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 976 ENDIF 977 978 IF( lwm ) THEN 979 ! write variables 980 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 981 ! close file 982 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 983 ENDIF 984 985 END SUBROUTINE rhg_cvg 821 986 822 987 … … 868 1033 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 869 1034 ! 870 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i 871 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i 1035 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) 1036 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 872 1037 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 873 1038 ! … … 876 1041 END SUBROUTINE rhg_evp_rst 877 1042 1043 878 1044 #else 879 1045 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
