- Timestamp:
- 2020-10-22T20:49:56+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11842_SI3-10_EAP
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_eap.F90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP
- Property svn:externals
-
old new 1 ^/utils/build/arch@HEAD arch 2 ^/utils/build/makenemo@HEAD makenemo 3 ^/utils/build/mk@HEAD mk 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 ^/vendors/FCM@HEAD ext/FCM 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 1 ^/utils/build/arch@12130 arch 2 ^/utils/build/makenemo@12191 makenemo 3 ^/utils/build/mk@11662 mk 4 ^/utils/tools_r4.0-HEAD@12672 tools 5 ^/vendors/AGRIF/dev@10586 ext/AGRIF 6 ^/vendors/FCM@10134 ext/FCM 7 ^/vendors/IOIPSL@9655 ext/IOIPSL 8 9 # SETTE mapping (inactive) 10 #^/utils/CI/sette@12135 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r13574 r13662 43 43 USE prtctl ! Print control 44 44 45 USE netcdf ! NetCDF library for convergence test 45 46 IMPLICIT NONE 46 47 PRIVATE … … 49 50 PUBLIC rhg_eap_rst ! called by icedyn_rhg.F90 50 51 51 REAL(wp), PARAMETER 52 REAL(wp), PARAMETER :: pphi = 3.141592653589793_wp/12._wp ! diamond shaped floe smaller angle (default phi = 30 deg) 52 53 53 54 ! look-up table for calculating structure tensor … … 60 61 !! * Substitutions 61 62 # include "vectopt_loop_substitute.h90" 63 64 !! for convergence tests 65 INTEGER :: ncvgid ! netcdf file id 66 INTEGER :: nvarid ! netcdf variable id 67 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 62 68 !!---------------------------------------------------------------------- 63 69 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 131 137 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 132 138 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 133 REAL(wp) :: zalph1, z1_alph1 ! alpha coef from Bouillon 2009 or Kimmritz 2017 139 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 140 REAl(wp) :: zbetau, zbetav 134 141 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 135 REAL(wp) :: zd elta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2, zdsT! temporary scalars142 REAL(wp) :: zds2, zdt, zdt2, zdiv, zdiv2, zdsT ! temporary scalars 136 143 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 137 144 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 138 145 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 139 146 ! 140 REAL(wp) :: zresm ! Maximal error on ice velocity141 147 REAL(wp) :: zintb, zintn ! dummy argument 142 148 REAL(wp) :: zfac_x, zfac_y 143 149 REAL(wp) :: zshear, zdum1, zdum2 144 REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components 145 REAL(wp) :: zalphar, zalphas ! for mechanical redistribution 146 REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution 147 ! 148 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 149 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 150 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 150 REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components 151 REAL(wp) :: zalphar, zalphas ! for mechanical redistribution 152 REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution 153 ! 154 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 155 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 156 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 157 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 151 158 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 152 159 ! … … 159 166 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 160 167 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 161 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence162 168 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 163 169 ! ! ocean surface (ssh_m) if ice is not embedded … … 173 179 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 174 180 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 175 REAL(wp), DIMENSION(jpi,jpj) :: zfmask , zwf! mask at F points for the ice181 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 176 182 177 183 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 178 184 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 179 185 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 186 !! --- check convergence 187 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 180 188 !! --- diags 181 REAL(wp) , DIMENSION(jpi,jpj) :: zmsk00182 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig 1, zsig2, zsig3189 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 190 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 183 191 !! --- SIMIP diags 184 192 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 191 199 192 200 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 201 ! 202 ! for diagnostics and convergence tests 203 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 207 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 208 END DO 209 END DO 193 210 ! 194 211 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... … … 205 222 206 223 ! Lateral boundary conditions on velocity (modify zfmask) 207 zwf(:,:) = zfmask(:,:)208 224 DO jj = 2, jpjm1 209 225 DO ji = fs_2, fs_jpim1 ! vector opt. 210 226 IF( zfmask(ji,jj) == 0._wp ) THEN 211 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) ) ) 227 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 228 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 212 229 ENDIF 213 230 END DO … … 215 232 DO jj = 2, jpjm1 216 233 IF( zfmask(1,jj) == 0._wp ) THEN 217 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )234 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 218 235 ENDIF 219 236 IF( zfmask(jpi,jj) == 0._wp ) THEN 220 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )237 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 221 238 ENDIF 222 239 END DO 223 240 DO ji = 2, jpim1 224 241 IF( zfmask(ji,1) == 0._wp ) THEN 225 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )242 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 226 243 ENDIF 227 244 IF( zfmask(ji,jpj) == 0._wp ) THEN 228 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )245 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 229 246 ENDIF 230 247 END DO … … 240 257 z1_ecc2 = 1._wp / ecc2 241 258 242 ! Time step for subcycling243 zdtevp = rdt_ice / REAL( nn_nevp )244 z1_dtevp = 1._wp / zdtevp245 246 259 ! alpha parameters (Bouillon 2009) 247 260 IF( .NOT. ln_aEVP ) THEN 248 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 261 zdtevp = rdt_ice / REAL( nn_nevp ) 262 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 263 zalph2 = zalph1 * z1_ecc2 264 249 265 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 250 ENDIF 266 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 267 ELSE 268 zdtevp = rdt_ice 269 ! zalpha parameters set later on adaptatively 270 ENDIF 271 z1_dtevp = 1._wp / zdtevp 251 272 252 273 ! Initialise stress tensor … … 264 285 265 286 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 266 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile287 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 267 288 ELSE ; zkt = 0._wp 268 289 ENDIF … … 335 356 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) 336 357 ! ice-bottom stress at U points 337 zvCr = zaU(ji,jj) * rn_ depfra * hu_n(ji,jj)338 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )358 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 359 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 339 360 ! ice-bottom stress at V points 340 zvCr = zaV(ji,jj) * rn_ depfra * hv_n(ji,jj)341 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )361 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 362 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 342 363 ! ice_bottom stress at T points 343 zvCr = at_i(ji,jj) * rn_ depfra * ht_n(ji,jj)344 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )364 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 365 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 345 366 END DO 346 367 END DO … … 365 386 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 366 387 ! 367 !!$ IF(ln_ctl) THEN ! Convergence test 368 !!$ DO jj = 1, jpjm1 369 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 370 !!$ zv_ice(:,jj) = v_ice(:,jj) 371 !!$ END DO 372 !!$ ENDIF 388 ! convergence test 389 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 390 DO jj = 1, jpj 391 DO ji = 1, jpi 392 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 393 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 394 END DO 395 END DO 396 ENDIF 373 397 374 398 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 375 DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points399 DO jj = 1, jpjm1 376 400 DO ji = 1, jpim1 377 401 … … 383 407 END DO 384 408 END DO 385 CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1. ) 386 387 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 388 DO ji = 2, jpi ! no vector loop 389 390 ! shear at T points 391 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 392 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 393 ! & ) * 0.25_wp * r1_e1e2t(ji,jj) .25 was an error 394 & ) * r1_e1e2t(ji,jj) 395 !WRITE(numout,*) 'shear output', ji, jj, zdsT 396 409 CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1._wp ) 410 411 DO jj = 2, jpjm1 412 DO ji = 2, jpim1 ! no vector loop 397 413 398 414 ! shear**2 at T points (doc eq. A16) … … 412 428 & ) * r1_e1e2t(ji,jj) 413 429 zdt2 = zdt * zdt 430 431 ! delta at T points 432 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 433 434 END DO 435 END DO 436 CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1._wp ) 414 437 438 ! P/delta at T points 439 DO jj = 1, jpj 440 DO ji = 1, jpi 441 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 442 END DO 443 END DO 444 445 DO jj = 2, jpj ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 446 DO ji = 2, jpi ! no vector loop 447 448 ! shear at T points 449 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 450 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 451 & ) * 0.25_wp * r1_e1e2t(ji,jj) 452 !WRITE(numout,*) 'shear output', ji, jj, zdsT 453 454 ! divergence at T points (duplication to avoid communications) 455 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 456 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 457 & ) * r1_e1e2t(ji,jj) 458 459 ! tension at T points (duplication to avoid communications) 460 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) & 461 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 462 & ) * r1_e1e2t(ji,jj) 463 415 464 ! --- anisotropic stress calculation --- ! 416 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 417 zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 418 zstressptmp, zstressmtmp, & 419 zstress12tmp(ji,jj), strength(ji,jj), & 420 zalphar, zalphas) 465 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 466 zstressptmp, zstressmtmp, zstress12tmp(ji,jj), strength(ji,jj), zalphar, zalphas) 421 467 422 468 ! structure tensor update 423 IF (mod(jter,10) == 0) THEN 424 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 425 paniso_11(ji,jj), paniso_12(ji,jj), & 426 zmresult11, zmresult12) 427 428 paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth - zmresult11) * z1dtevpkth ! implicit 429 paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A - zmresult12) * z1dtevpkth ! implicit 430 ENDIF 431 469 !!$ IF (mod(jter,10) == 0) THEN 470 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11, zmresult12) 471 472 !!$ paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth + zmresult11) * z1dtevpkth ! implicit 473 !!$ paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A + zmresult12) * z1dtevpkth ! implicit 474 paniso_11(ji,jj) = (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 475 paniso_12(ji,jj) = (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 476 !!$ ENDIF 432 477 433 478 IF (jter == nn_nevp) THEN … … 438 483 ENDIF 439 484 440 ! delta at T points 441 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 442 443 ! P/delta at T points 444 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 445 446 ! alpha & beta for aEVP 485 ! alpha for aEVP 447 486 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 448 487 ! alpha = beta = sqrt(4*gamma) 449 !IF( ln_aEVP ) THEN 450 ! zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 451 ! z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 452 !ENDIF 453 488 IF( ln_aEVP ) THEN 489 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 490 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 491 zalph2 = zalph1 492 z1_alph2 = z1_alph1 493 ! explicit: 494 ! z1_alph1 = 1._wp / zalph1 495 ! z1_alph2 = 1._wp / zalph1 496 ! zalph1 = zalph1 - 1._wp 497 ! zalph2 = zalph1 498 ENDIF 499 454 500 ! stress at T points (zkt/=0 if landfast) 455 !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 456 !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 457 zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 458 zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 501 !zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 502 !zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 503 !!$ zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 504 !!$ zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 505 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1 506 zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1 459 507 !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1 ) * z1_alph1 460 508 !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1 ) * z1_alph1 … … 462 510 END DO 463 511 END DO 464 !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )465 512 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 513 514 ! Save beta at T-points for further computations 515 IF( ln_aEVP ) THEN 516 DO jj = 1, jpj 517 DO ji = 1, jpi 518 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 519 END DO 520 END DO 521 ENDIF 466 522 467 523 DO jj = 1, jpjm1 … … 469 525 ! stress12tmp at F points 470 526 zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 471 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) &472 & ) * 0.25_wp * r1_e1e2f(ji,jj)473 474 ! alpha & betafor aEVP527 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & 528 & ) * 0.25_wp * r1_e1e2f(ji,jj) 529 530 ! alpha for aEVP 475 531 IF( ln_aEVP ) THEN 476 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 477 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 532 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 533 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 534 ! explicit: 535 ! z1_alph2 = 1._wp / zalph2 536 ! zalph2 = zalph2 - 1._wp 478 537 ENDIF 479 538 480 539 ! stress at F points (zkt/=0 if landfast) 481 540 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 482 zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 541 !!$ zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 542 zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 483 543 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1 484 544 485 545 END DO 486 546 END DO 487 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )547 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 488 548 489 549 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! … … 545 605 ! 546 606 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 547 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 548 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 549 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 550 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 551 & ) * 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 607 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 608 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 609 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 610 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 611 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 612 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 613 & ) / ( zbetav + 1._wp ) & 614 & ) * 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 552 615 & ) * zmsk00y(ji,jj) 553 616 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 554 v_ice(ji,jj) = ( ( 555 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)556 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast557 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0558 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) 617 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 618 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 619 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 620 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 621 & ) * 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 559 622 & ) * zmsk00y(ji,jj) 560 623 ENDIF … … 596 659 ! 597 660 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 598 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 599 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 600 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 601 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 602 & ) * 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 661 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 662 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 663 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 664 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 665 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 666 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 667 & ) / ( zbetau + 1._wp ) & 668 & ) * 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 603 669 & ) * zmsk00x(ji,jj) 604 670 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 605 u_ice(ji,jj) = ( ( 606 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)607 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast608 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0609 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) 671 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 672 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 673 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 674 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 675 & ) * 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 610 676 & ) * zmsk00x(ji,jj) 611 677 ENDIF … … 649 715 ! 650 716 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 651 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 652 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 653 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 654 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 655 & ) * 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 717 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 718 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 719 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 720 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 721 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 722 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 723 & ) / ( zbetau + 1._wp ) & 724 & ) * 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 656 725 & ) * zmsk00x(ji,jj) 657 726 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 658 u_ice(ji,jj) = ( ( 659 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)660 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast661 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0662 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) 727 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 728 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 729 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 730 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 731 & ) * 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 663 732 & ) * zmsk00x(ji,jj) 664 733 ENDIF … … 700 769 ! 701 770 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 702 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 703 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 704 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 705 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 706 & ) * 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 771 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 772 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 773 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 774 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 775 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 776 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 777 & ) / ( zbetav + 1._wp ) & 778 & ) * 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 707 779 & ) * zmsk00y(ji,jj) 708 780 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 709 v_ice(ji,jj) = ( ( 710 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)711 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast712 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0713 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) 781 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 782 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 783 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 784 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 785 & ) * 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 714 786 & ) * zmsk00y(ji,jj) 715 787 ENDIF … … 726 798 ENDIF 727 799 728 !!$ IF(ln_ctl) THEN ! Convergence test 729 !!$ DO jj = 2 , jpjm1 730 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 731 !!$ END DO 732 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 733 !!$ CALL mpp_max( 'icedyn_rhg_eap', zresm ) ! max over the global domain 734 !!$ ENDIF 800 ! convergence test 801 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 735 802 ! 736 737 803 ! ! ==================== ! 738 804 END DO ! end loop over jter ! 739 805 ! ! ==================== ! 806 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 740 807 ! 741 808 CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1. ) ! only need this in ridging module after rheology completed … … 764 831 zdt2 = zdt * zdt 765 832 833 zten_i(ji,jj) = zdt 834 766 835 ! shear**2 at T points (doc eq. A16) 767 836 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 778 847 779 848 ! delta at T points 780 z delta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )781 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0782 pdelta_i(ji,jj) = z delta + rn_creepl * rswitch849 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 850 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 851 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 783 852 784 853 END DO 785 854 END DO 786 CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )787 !CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1. )855 CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 856 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. ) 788 857 789 858 ! --- Store the stress tensor for the next time step --- ! 790 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )791 859 pstress1_i (:,:) = zs1 (:,:) 792 860 pstress2_i (:,:) = zs2 (:,:) … … 797 865 ! 5) diagnostics 798 866 !------------------------------------------------------------------------------! 799 DO jj = 1, jpj800 DO ji = 1, jpi801 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice802 END DO803 END DO804 805 867 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 806 868 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 821 883 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 822 884 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 885 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 823 886 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 824 887 825 ! --- stress tensor--- !826 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN888 ! --- Stress tensor invariants (SIMIP diags) --- ! 889 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 827 890 ! 828 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )891 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 829 892 ! 830 DO jj = 2, jpjm1 831 DO ji = 2, jpim1 832 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 833 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 834 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 835 836 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 837 838 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 839 840 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 841 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 842 !! 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 843 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 844 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 845 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 846 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 893 DO jj = 1, jpj 894 DO ji = 1, jpi 895 896 ! Ice stresses 897 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 898 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 899 ! I know, this can be confusing... 900 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 901 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 902 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 903 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 904 905 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 906 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 907 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 908 847 909 END DO 848 END DO 849 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 910 END DO 850 911 ! 851 CALL iom_put( 'isig1' , zsig1 ) 852 CALL iom_put( 'isig2' , zsig2 ) 853 CALL iom_put( 'isig3' , zsig3 ) 912 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 913 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 914 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 915 916 DEALLOCATE ( zsig_I, zsig_II ) 917 918 ENDIF 919 920 ! --- Normalized stress tensor principal components --- ! 921 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 922 ! Recommendation 1 : we use ice strength, not replacement pressure 923 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 924 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 854 925 ! 855 ! Stress tensor invariants (normal and shear stress N/m) 856 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 857 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 858 859 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 926 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 927 ! 928 DO jj = 1, jpj 929 DO ji = 1, jpi 930 931 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 932 ! and **deformations** at current iterates 933 ! following Lemieux & Dupont (2020) 934 zfac = zp_delt(ji,jj) 935 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 936 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 937 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 938 939 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 940 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 941 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 942 943 ! Normalized principal stresses (used to display the ellipse) 944 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 945 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 946 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 947 END DO 948 END DO 949 ! 950 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 951 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 952 953 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 954 860 955 ENDIF 861 956 … … 931 1026 ENDIF 932 1027 ! 1028 ! --- convergence tests --- ! 1029 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 1030 IF( iom_use('uice_cvg') ) THEN 1031 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 1032 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 1033 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 1034 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 1035 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 1036 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 1037 ENDIF 1038 ENDIF 1039 ENDIF 1040 ! 1041 DEALLOCATE( zmsk00, zmsk15 ) 1042 ! 933 1043 END SUBROUTINE ice_dyn_rhg_eap 934 1044 935 SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 936 pshear, pa11, pa12, & 937 pstressp, pstressm, & 938 pstress12, strength, & 939 palphar, palphas) 1045 1046 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 1047 !!---------------------------------------------------------------------- 1048 !! *** ROUTINE rhg_cvg *** 1049 !! 1050 !! ** Purpose : check convergence of oce rheology 1051 !! 1052 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 1053 !! during the sub timestepping of rheology so as: 1054 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 1055 !! This routine is called every sub-iteration, so it is cpu expensive 1056 !! 1057 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1058 !!---------------------------------------------------------------------- 1059 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1060 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 1061 !! 1062 INTEGER :: it, idtime, istatus 1063 INTEGER :: ji, jj ! dummy loop indices 1064 REAL(wp) :: zresm ! local real 1065 CHARACTER(len=20) :: clname 1066 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 1067 !!---------------------------------------------------------------------- 1068 1069 ! create file 1070 IF( kt == nit000 .AND. kiter == 1 ) THEN 1071 ! 1072 IF( lwp ) THEN 1073 WRITE(numout,*) 1074 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 1075 WRITE(numout,*) '~~~~~~~' 1076 ENDIF 1077 ! 1078 IF( lwm ) THEN 1079 clname = 'ice_cvg.nc' 1080 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 1081 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 1082 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 1083 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 1084 istatus = NF90_ENDDEF(ncvgid) 1085 ENDIF 1086 ! 1087 ENDIF 1088 1089 ! time 1090 it = ( kt - 1 ) * kitermax + kiter 1091 1092 ! convergence 1093 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 1094 zresm = 0._wp 1095 ELSE 1096 DO jj = 1, jpj 1097 DO ji = 1, jpi 1098 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1099 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 1100 END DO 1101 END DO 1102 zresm = MAXVAL( zres ) 1103 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1104 ENDIF 1105 1106 IF( lwm ) THEN 1107 ! write variables 1108 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1109 ! close file 1110 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 1111 ENDIF 1112 1113 END SUBROUTINE rhg_cvg 1114 1115 1116 SUBROUTINE update_stress_rdg( ksub, kndte, pdivu, ptension, pshear, pa11, pa12, & 1117 & pstressp, pstressm, pstress12, pstrength, palphar, palphas ) 940 1118 !!--------------------------------------------------------------------- 941 1119 !! *** SUBROUTINE update_stress_rdg *** … … 945 1123 !!--------------------------------------------------------------------- 946 1124 INTEGER, INTENT(in ) :: ksub, kndte 947 REAL(wp), INTENT(in ) :: strength1125 REAL(wp), INTENT(in ) :: pstrength 948 1126 REAL(wp), INTENT(in ) :: pdivu, ptension, pshear 949 1127 REAL(wp), INTENT(in ) :: pa11, pa12 … … 992 1170 IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN 993 1171 zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 994 4._wp*pa12*pa12 )1172 4._wp*pa12*pa12 ) 995 1173 996 1174 zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma !Cos^2 … … 1067 1245 ! % range in kx 1068 1246 kx = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 1069 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1247 !!clem kx = MAX( 1, MIN( nx_yield-1, INT((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 ) ) 1248 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1070 1249 1071 1250 ky = int(zy*zinvdy) + 1 1251 !!clem ky = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) ) 1072 1252 kyw = ky - zy*zinvdy 1073 1253 1074 1254 ka = int((zatempprime-0.5_wp)*zinvda) + 1 1255 !!clem ka = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) ) 1075 1256 kaw = ka - (zatempprime-0.5_wp)*zinvda 1076 1257 1077 1258 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 1078 zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1079 & + (1._wp-zkxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1080 & + zkxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1081 & + zkxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1082 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1083 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1084 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1085 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1086 zstemp12r = zkxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1087 & + (1._wp-zkxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1088 & + zkxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1089 & + zkxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1090 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1091 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1092 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1093 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1094 zstemp22r = zkxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1095 & + (1._wp-zkxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1096 & + zkxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1097 & + zkxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1098 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1099 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1100 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1101 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1259 !!$ zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1260 !!$ & + (1._wp-zkxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1261 !!$ & + zkxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1262 !!$ & + zkxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1263 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1264 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1265 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1266 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1267 !!$ zstemp12r = zkxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1268 !!$ & + (1._wp-zkxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1269 !!$ & + zkxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1270 !!$ & + zkxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1271 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1272 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1273 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1274 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1275 !!$ zstemp22r = zkxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1276 !!$ & + (1._wp-zkxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1277 !!$ & + zkxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1278 !!$ & + zkxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1279 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1280 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1281 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1282 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1283 !!$ 1284 !!$ zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1285 !!$ & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1286 !!$ & + zkxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1287 !!$ & + zkxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1288 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1289 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1290 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1291 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1292 !!$ zstemp12s = zkxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1293 !!$ & + (1._wp-zkxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1294 !!$ & + zkxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1295 !!$ & + zkxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1296 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1297 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1298 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1299 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1300 !!$ zstemp22s = zkxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1301 !!$ & + (1._wp-zkxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1302 !!$ & + zkxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1303 !!$ & + zkxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1304 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1305 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1306 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1307 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1308 1309 zstemp11r = s11r(kx,ky,ka) 1310 zstemp12r = s12r(kx,ky,ka) 1311 zstemp22r = s22r(kx,ky,ka) 1312 zstemp11s = s11s(kx,ky,ka) 1313 zstemp12s = s12s(kx,ky,ka) 1314 zstemp22s = s22s(kx,ky,ka) 1102 1315 1103 zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1104 & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1105 & + zkxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1106 & + zkxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1107 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1108 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1109 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1110 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1111 zstemp12s = zkxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1112 & + (1._wp-zkxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1113 & + zkxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1114 & + zkxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1115 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1116 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1117 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1118 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1119 zstemp22s = zkxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1120 & + (1._wp-zkxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1121 & + zkxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1122 & + zkxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1123 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1124 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1125 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1126 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1127 1316 1128 1317 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1129 1318 ! Tsamados 2013) 1130 1319 1131 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin1132 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin1133 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin1320 zsig11 = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin 1321 zsig12 = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin 1322 zsig22 = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1134 1323 1135 1324 !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 … … 1138 1327 1139 1328 ! Update stress 1140 zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 - 1329 zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 - 2._wp*zQ11Q12 *zsig12 1141 1330 zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12 1142 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 1331 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 2._wp*zQ11Q12 *zsig12 1143 1332 1144 1333 pstressp = zsgprm11 + zsgprm22 … … 1150 1339 IF (ksub == kndte) THEN 1151 1340 zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 1152 + zQ12Q12*zstemp22r1153 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) &1154 - zQ12Q12*zstemp12r1341 + zQ12Q12*zstemp22r 1342 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) & 1343 - zQ12Q12*zstemp12r 1155 1344 zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 1156 + zQ11Q11*zstemp22r1345 + zQ11Q11*zstemp22r 1157 1346 1158 1347 zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s & 1159 + zQ12Q12*zstemp22s1160 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) &1161 - zQ12Q12*zstemp12s1348 + zQ12Q12*zstemp22s 1349 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) & 1350 - zQ12Q12*zstemp12s 1162 1351 zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 1163 + zQ11Q11*zstemp22s1352 + zQ11Q11*zstemp22s 1164 1353 1165 1354 palphar = zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 1166 + zrotstemp22r*zdtemp221355 + zrotstemp22r*zdtemp22 1167 1356 palphas = zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 1168 + zrotstemp22s*zdtemp221357 + zrotstemp22s*zdtemp22 1169 1358 ENDIF 1170 1359 END SUBROUTINE update_stress_rdg … … 1173 1362 1174 1363 1175 SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 1176 pa11, pa12, & 1177 pmresult11, pmresult12) 1364 SUBROUTINE calc_ffrac( pstressp, pstressm, pstress12, pa11, pa12, & 1365 & pmresult11, pmresult12 ) 1178 1366 !!--------------------------------------------------------------------- 1179 1367 !! *** ROUTINE calc_ffrac *** … … 1193 1381 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1194 1382 1195 REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation 1383 !!$ REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation 1384 REAL (wp), PARAMETER :: kfrac = 1.e-3_wp ! rate of fracture formation 1196 1385 REAL (wp), PARAMETER :: threshold = 0.3_wp ! critical confinement ratio 1197 1386 !!--------------------------------------------------------------- … … 1224 1413 ENDIF 1225 1414 1226 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 & 1227 + zQ12Q12*zsigma22 ! S(1,1) 1228 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 & 1229 + zQ11Q11*zsigma22 ! S(2,2) 1415 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 + zQ12Q12*zsigma22 ! S(1,1) 1416 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 + zQ11Q11*zsigma22 ! S(2,2) 1230 1417 1231 1418 ! Pure divergence … … 1238 1425 ! which leads to the loss of their shape, so we again model it through diffusion 1239 1426 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1240 pmresult11 = kfrac * (pa11 - zQ12Q12)1241 pmresult12 = kfrac * (pa12 + zQ11Q12)1427 pmresult11 = - kfrac * (pa11 - zQ12Q12) 1428 pmresult12 = - kfrac * (pa12 + zQ11Q12) 1242 1429 1243 1430 ! Shear faulting … … 1246 1433 pmresult12 = 0.0_wp 1247 1434 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 1248 pmresult11 = kfrac * (pa11 - zQ12Q12)1249 pmresult12 = kfrac * (pa12 + zQ11Q12)1435 pmresult11 = - kfrac * (pa11 - zQ12Q12) 1436 pmresult12 = - kfrac * (pa12 + zQ11Q12) 1250 1437 1251 1438 ! Horizontal spalling … … 1278 1465 REAL(wp) :: da, dx, dy, dp, dz, a1 1279 1466 1467 !!clem 1468 REAL(wp) :: zw1, zw2, zfac, ztemp 1469 REAL(wp) :: idx, idy, idz 1470 1280 1471 REAL(wp), PARAMETER :: eps6 = 1.0e-6_wp 1281 1472 !!---------------------------------------------------------------------- … … 1288 1479 id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) 1289 1480 id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) 1290 id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. )1291 id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. )1481 id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. ) 1482 id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. ) 1292 1483 ! 1293 1484 IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN ! fields exist … … 1295 1486 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) 1296 1487 CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 1297 CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11)1298 CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12)1488 CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11 ) 1489 CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12 ) 1299 1490 ELSE ! start rheology from rest 1300 1491 IF(lwp) WRITE(numout,*) … … 1303 1494 stress2_i (:,:) = 0._wp 1304 1495 stress12_i(:,:) = 0._wp 1305 aniso_11 (:,:) = 0.5_wp1306 aniso_12 (:,:) = 0._wp1496 aniso_11 (:,:) = 0.5_wp 1497 aniso_12 (:,:) = 0._wp 1307 1498 ENDIF 1308 1499 ELSE !* Start from rest … … 1312 1503 stress2_i (:,:) = 0._wp 1313 1504 stress12_i(:,:) = 0._wp 1314 aniso_11 (:,:) = 0.5_wp1315 aniso_12 (:,:) = 0._wp1505 aniso_11 (:,:) = 0.5_wp 1506 aniso_12 (:,:) = 0._wp 1316 1507 ENDIF 1317 1508 ! 1318 1509 1319 da = 0.5_wp/real(na_yield-1,kind=wp) 1320 ainit = 0.5_wp - da 1321 dx = rpi/real(nx_yield-1,kind=wp) 1322 xinit = rpi + 0.25_wp*rpi - dx 1323 dz = rpi/real(nz,kind=wp) 1324 zinit = -rpi*0.5_wp 1325 dy = rpi/real(ny_yield-1,kind=wp) 1326 yinit = -dy 1327 1328 DO ia=1,na_yield 1329 DO ix=1,nx_yield 1330 DO iy=1,ny_yield 1331 s11r(ix,iy,ia) = 0._wp 1332 s12r(ix,iy,ia) = 0._wp 1333 s22r(ix,iy,ia) = 0._wp 1334 s11s(ix,iy,ia) = 0._wp 1335 s12s(ix,iy,ia) = 0._wp 1336 s22s(ix,iy,ia) = 0._wp 1337 IF (ia <= na_yield-1) THEN 1338 DO iz=1,nz 1339 s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1340 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1341 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1342 s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1343 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1344 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1345 s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1346 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1347 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1348 s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1349 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1350 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1351 s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1352 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1353 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1354 s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1355 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1356 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1357 ENDDO 1358 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1359 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1360 IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1361 IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1362 IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1363 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1364 ELSE 1365 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1366 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1367 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1368 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1369 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1370 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1371 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1372 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1373 IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1374 IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1375 IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1376 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1377 ENDIF 1378 ENDDO 1379 ENDDO 1380 ENDDO 1510 da = 0.5_wp/real(na_yield-1,kind=wp) 1511 ainit = 0.5_wp - da 1512 dx = rpi/real(nx_yield-1,kind=wp) 1513 xinit = rpi + 0.25_wp*rpi - dx 1514 dz = rpi/real(nz,kind=wp) 1515 zinit = -rpi*0.5_wp 1516 dy = rpi/real(ny_yield-1,kind=wp) 1517 yinit = -dy 1518 1519 s11r(:,:,:) = 0._wp 1520 s12r(:,:,:) = 0._wp 1521 s22r(:,:,:) = 0._wp 1522 s11s(:,:,:) = 0._wp 1523 s12s(:,:,:) = 0._wp 1524 s22s(:,:,:) = 0._wp 1525 1526 !!$ DO ia=1,na_yield 1527 !!$ DO ix=1,nx_yield 1528 !!$ DO iy=1,ny_yield 1529 !!$ s11r(ix,iy,ia) = 0._wp 1530 !!$ s12r(ix,iy,ia) = 0._wp 1531 !!$ s22r(ix,iy,ia) = 0._wp 1532 !!$ s11s(ix,iy,ia) = 0._wp 1533 !!$ s12s(ix,iy,ia) = 0._wp 1534 !!$ s22s(ix,iy,ia) = 0._wp 1535 !!$ IF (ia <= na_yield-1) THEN 1536 !!$ DO iz=1,nz 1537 !!$ s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1538 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1539 !!$ s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1540 !!$ s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1541 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1542 !!$ s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1543 !!$ s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1544 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1545 !!$ s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1546 !!$ s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1547 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1548 !!$ s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1549 !!$ s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1550 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1551 !!$ s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1552 !!$ s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1553 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1554 !!$ s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1555 !!$ ENDDO 1556 !!$ IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1557 !!$ IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1558 !!$ IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1559 !!$ IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1560 !!$ IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1561 !!$ IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1562 !!$ ELSE 1563 !!$ s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1564 !!$ s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1565 !!$ s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1566 !!$ s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1567 !!$ s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1568 !!$ s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1569 !!$ IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1570 !!$ IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1571 !!$ IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1572 !!$ IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1573 !!$ IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1574 !!$ IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1575 !!$ ENDIF 1576 !!$ ENDDO 1577 !!$ ENDDO 1578 !!$ ENDDO 1579 1580 !! faster but still very slow => to be improved 1581 zfac = dz/sin(2._wp*pphi) 1582 DO ia = 1, na_yield-1 1583 zw1 = w1(ainit+ia*da) 1584 zw2 = w2(ainit+ia*da) 1585 DO iz = 1, nz 1586 idz = zinit+iz*dz 1587 ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 1588 DO iy = 1, ny_yield 1589 idy = yinit+iy*dy 1590 DO ix = 1, nx_yield 1591 idx = xinit+ix*dx 1592 s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac 1593 s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac 1594 s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac 1595 s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac 1596 s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac 1597 s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac 1598 END DO 1599 END DO 1600 END DO 1601 END DO 1602 1603 zfac = 1._wp/sin(2._wp*pphi) 1604 ia = na_yield 1605 DO iy = 1, ny_yield 1606 idy = yinit+iy*dy 1607 DO ix = 1, nx_yield 1608 idx = xinit+ix*dx 1609 s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac 1610 s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac 1611 s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac 1612 s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac 1613 s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac 1614 s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac 1615 ENDDO 1616 ENDDO 1617 WHERE (ABS(s11r(:,:,:)) < eps6) s11r(:,:,:) = 0._wp 1618 WHERE (ABS(s12r(:,:,:)) < eps6) s12r(:,:,:) = 0._wp 1619 WHERE (ABS(s22r(:,:,:)) < eps6) s22r(:,:,:) = 0._wp 1620 WHERE (ABS(s11s(:,:,:)) < eps6) s11s(:,:,:) = 0._wp 1621 WHERE (ABS(s12s(:,:,:)) < eps6) s12s(:,:,:) = 0._wp 1622 WHERE (ABS(s22s(:,:,:)) < eps6) s22s(:,:,:) = 0._wp 1623 1381 1624 1382 1625 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file … … 1388 1631 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 1389 1632 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 1390 CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11)1391 CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12)1633 CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11 ) 1634 CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12 ) 1392 1635 ! 1393 1636 ENDIF … … 1404 1647 !!------------------------------------------------------------------- 1405 1648 1406 w1 = - 223.87569446_wp &1407 & + 2361.2198663_wp*pa &1649 w1 = - 223.87569446_wp & 1650 & + 2361.21986630_wp*pa & 1408 1651 & - 10606.56079975_wp*pa*pa & 1409 1652 & + 26315.50025642_wp*pa*pa*pa & … … 1411 1654 & + 34397.72407466_wp*pa*pa*pa*pa*pa & 1412 1655 & - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 1413 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa1656 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 1414 1657 1415 1658 END FUNCTION w1 … … 1424 1667 !!------------------------------------------------------------------- 1425 1668 1426 w2 = - 6670.68911883_wp &1427 & + 70222.33061536_wp*pa &1428 & - 314871.71525448_wp*pa*pa &1429 & + 779570.02793492_wp*pa*pa*pa &1669 w2 = - 6670.68911883_wp & 1670 & + 70222.33061536_wp*pa & 1671 & - 314871.71525448_wp*pa*pa & 1672 & + 779570.02793492_wp*pa*pa*pa & 1430 1673 & - 1151098.82436864_wp*pa*pa*pa*pa & 1431 1674 & + 1013896.59464498_wp*pa*pa*pa*pa*pa & 1432 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa &1433 & + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa1675 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 1676 & + 102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa 1434 1677 1435 1678 END FUNCTION w2 … … 1789 2032 !! Default option Empty module NO SI3 sea-ice model 1790 2033 !!---------------------------------------------------------------------- 1791 1792 2034 #endif 2035 1793 2036 !!============================================================================== 1794 2037 END MODULE icedyn_rhg_eap
Note: See TracChangeset
for help on using the changeset viewer.
