Changeset 13746 for NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
- Timestamp:
- 2020-11-08T22:25:28+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
r13574 r13746 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 REAL(wp) :: invw ! for test case 148 ! 149 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 150 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 151 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 REAL(wp) :: zinvw ! for test case 154 155 ! 156 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 157 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 158 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 159 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 152 160 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 153 161 ! … … 160 168 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 161 169 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 162 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence163 170 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 164 171 ! ! ocean surface (ssh_m) if ice is not embedded … … 174 181 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 175 182 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 176 REAL(wp), DIMENSION(jpi,jpj) :: zfmask , zwf! mask at F points for the ice183 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 177 184 178 185 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 179 186 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 180 187 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 188 !! --- check convergence 189 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 181 190 !! --- diags 182 REAL(wp) , DIMENSION(jpi,jpj) :: zmsk00183 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig 1, zsig2, zsig3191 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 192 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 184 193 !! --- SIMIP diags 185 194 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 192 201 193 202 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 203 ! 204 ! for diagnostics and convergence tests 205 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 206 DO jj = 1, jpj 207 DO ji = 1, jpi 208 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 209 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 210 END DO 211 END DO 194 212 ! 195 213 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... … … 206 224 207 225 ! Lateral boundary conditions on velocity (modify zfmask) 208 zwf(:,:) = zfmask(:,:)209 226 DO jj = 2, jpjm1 210 227 DO ji = fs_2, fs_jpim1 ! vector opt. 211 228 IF( zfmask(ji,jj) == 0._wp ) THEN 212 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) ) ) 229 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 230 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 213 231 ENDIF 214 232 END DO … … 216 234 DO jj = 2, jpjm1 217 235 IF( zfmask(1,jj) == 0._wp ) THEN 218 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )236 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 219 237 ENDIF 220 238 IF( zfmask(jpi,jj) == 0._wp ) THEN 221 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )239 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 222 240 ENDIF 223 241 END DO 224 242 DO ji = 2, jpim1 225 243 IF( zfmask(ji,1) == 0._wp ) THEN 226 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )244 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 227 245 ENDIF 228 246 IF( zfmask(ji,jpj) == 0._wp ) THEN 229 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )247 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 230 248 ENDIF 231 249 END DO … … 237 255 zrhoco = rau0 * rn_cio 238 256 !extra code for test case boundary conditions 239 invw=1._wp/(zrhoco*0.5_wp)257 zinvw=1._wp/(zrhoco*0.5_wp) 240 258 241 259 ! ecc2: square of yield ellipse eccenticrity … … 243 261 z1_ecc2 = 1._wp / ecc2 244 262 245 ! Time step for subcycling246 zdtevp = rdt_ice / REAL( nn_nevp )247 z1_dtevp = 1._wp / zdtevp248 249 263 ! alpha parameters (Bouillon 2009) 250 264 IF( .NOT. ln_aEVP ) THEN 251 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 265 zdtevp = rdt_ice / REAL( nn_nevp ) 266 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 267 zalph2 = zalph1 * z1_ecc2 268 252 269 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 253 ENDIF 270 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 271 ELSE 272 zdtevp = rdt_ice 273 ! zalpha parameters set later on adaptatively 274 ENDIF 275 z1_dtevp = 1._wp / zdtevp 254 276 255 277 ! Initialise stress tensor … … 267 289 268 290 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 269 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile291 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 270 292 ELSE ; zkt = 0._wp 271 293 ENDIF … … 338 360 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) 339 361 ! ice-bottom stress at U points 340 zvCr = zaU(ji,jj) * rn_ depfra * hu_n(ji,jj)341 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )362 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 363 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 342 364 ! ice-bottom stress at V points 343 zvCr = zaV(ji,jj) * rn_ depfra * hv_n(ji,jj)344 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )365 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 366 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 345 367 ! ice_bottom stress at T points 346 zvCr = at_i(ji,jj) * rn_ depfra * ht_n(ji,jj)347 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )368 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 369 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 348 370 END DO 349 371 END DO … … 368 390 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 369 391 ! 370 !!$ IF(ln_ctl) THEN ! Convergence test 371 !!$ DO jj = 1, jpjm1 372 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 373 !!$ zv_ice(:,jj) = v_ice(:,jj) 374 !!$ END DO 375 !!$ ENDIF 392 ! convergence test 393 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 394 DO jj = 1, jpj 395 DO ji = 1, jpi 396 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 397 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 398 END DO 399 END DO 400 ENDIF 376 401 377 402 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 378 DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points403 DO jj = 1, jpjm1 379 404 DO ji = 1, jpim1 380 405 … … 386 411 END DO 387 412 END DO 388 CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1. ) 389 390 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 391 DO ji = 2, jpi ! no vector loop 392 393 ! shear at T points 394 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 395 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 396 & ) * 0.25_wp * r1_e1e2t(ji,jj) 397 !WRITE(numout,*) 'shear output', ji, jj, zdsT 398 413 CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1._wp ) 414 415 DO jj = 2, jpjm1 416 DO ji = 2, jpim1 ! no vector loop 399 417 400 418 ! shear**2 at T points (doc eq. A16) … … 414 432 & ) * r1_e1e2t(ji,jj) 415 433 zdt2 = zdt * zdt 434 435 ! delta at T points 436 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 437 438 END DO 439 END DO 440 CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1._wp ) 416 441 442 ! P/delta at T points 443 DO jj = 1, jpj 444 DO ji = 1, jpi 445 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 446 END DO 447 END DO 448 449 DO jj = 2, jpj ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 450 DO ji = 2, jpi ! no vector loop 451 452 ! shear at T points 453 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 454 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 455 & ) * 0.25_wp * r1_e1e2t(ji,jj) 456 !WRITE(numout,*) 'shear output', ji, jj, zdsT 457 458 ! divergence at T points (duplication to avoid communications) 459 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 460 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 461 & ) * r1_e1e2t(ji,jj) 462 463 ! tension at T points (duplication to avoid communications) 464 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) & 465 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 466 & ) * r1_e1e2t(ji,jj) 467 417 468 ! --- anisotropic stress calculation --- ! 418 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 419 zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 420 zstressptmp, zstressmtmp, & 421 zstress12tmp(ji,jj), strength(ji,jj), & 422 zalphar, zalphas) 469 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 470 zstressptmp, zstressmtmp, zstress12tmp(ji,jj), strength(ji,jj), zalphar, zalphas) 423 471 424 472 ! structure tensor update 425 IF (mod(jter,10) == 0) THEN 426 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 427 paniso_11(ji,jj), paniso_12(ji,jj), & 428 zmresult11, zmresult12) 429 430 paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth - zmresult11) * z1dtevpkth ! implicit 431 paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A - zmresult12) * z1dtevpkth ! implicit 432 ENDIF 433 473 !!$ IF (mod(jter,10) == 0) THEN 474 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11, zmresult12) 475 476 !!$ paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth + zmresult11) * z1dtevpkth ! implicit 477 !!$ paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A + zmresult12) * z1dtevpkth ! implicit 478 paniso_11(ji,jj) = (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 479 paniso_12(ji,jj) = (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 480 !!$ ENDIF 434 481 435 482 IF (jter == nn_nevp) THEN … … 440 487 ENDIF 441 488 442 ! delta at T points 443 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 444 445 ! P/delta at T points 446 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 447 448 ! alpha & beta for aEVP 489 ! alpha for aEVP 449 490 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 450 491 ! alpha = beta = sqrt(4*gamma) 451 !IF( ln_aEVP ) THEN 452 ! zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 453 ! z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 454 !ENDIF 455 492 IF( ln_aEVP ) THEN 493 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 494 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 495 zalph2 = zalph1 496 z1_alph2 = z1_alph1 497 ! explicit: 498 ! z1_alph1 = 1._wp / zalph1 499 ! z1_alph2 = 1._wp / zalph1 500 ! zalph1 = zalph1 - 1._wp 501 ! zalph2 = zalph1 502 ENDIF 503 456 504 ! stress at T points (zkt/=0 if landfast) 457 !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 458 !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 459 zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 460 zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 505 !zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 506 !zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 507 !!$ zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 508 !!$ zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 509 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1 510 zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1 461 511 !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1 ) * z1_alph1 462 512 !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1 ) * z1_alph1 … … 464 514 END DO 465 515 END DO 466 !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )467 516 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 517 518 ! Save beta at T-points for further computations 519 IF( ln_aEVP ) THEN 520 DO jj = 1, jpj 521 DO ji = 1, jpi 522 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 523 END DO 524 END DO 525 ENDIF 468 526 469 527 DO jj = 1, jpjm1 … … 471 529 ! stress12tmp at F points 472 530 zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 473 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) &474 & ) * 0.25_wp * r1_e1e2f(ji,jj)475 476 ! alpha & betafor aEVP531 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & 532 & ) * 0.25_wp * r1_e1e2f(ji,jj) 533 534 ! alpha for aEVP 477 535 IF( ln_aEVP ) THEN 478 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 479 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 536 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 537 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 538 ! explicit: 539 ! z1_alph2 = 1._wp / zalph2 540 ! zalph2 = zalph2 - 1._wp 480 541 ENDIF 481 542 482 543 ! stress at F points (zkt/=0 if landfast) 483 544 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 484 zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 545 !!$ zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 546 zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 485 547 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1 486 548 487 549 END DO 488 550 END DO 489 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )551 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 490 552 491 553 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! … … 547 609 ! 548 610 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 549 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 550 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 551 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 552 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 553 & ) * 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 611 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 612 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 613 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 614 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 615 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 616 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 617 & ) / ( zbetav + 1._wp ) & 618 & ) * 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 554 619 & ) * zmsk00y(ji,jj) 555 620 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 556 v_ice(ji,jj) = ( ( 557 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)558 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast559 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0560 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) 621 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 622 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 623 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 624 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 625 & ) * 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 561 626 & ) * zmsk00y(ji,jj) 562 627 ENDIF 563 628 !extra code for test case boundary conditions 564 629 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 565 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))630 v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 566 631 END IF 632 567 633 END DO 568 634 END DO … … 602 668 ! 603 669 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 604 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 605 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 606 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 607 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 608 & ) * 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 670 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 671 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(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) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 674 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 675 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 676 & ) / ( zbetau + 1._wp ) & 677 & ) * 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 609 678 & ) * zmsk00x(ji,jj) 610 679 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 611 u_ice(ji,jj) = ( ( 612 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)613 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast614 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0615 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) 680 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 681 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 682 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 683 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 684 & ) * 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 616 685 & ) * zmsk00x(ji,jj) 617 686 ENDIF 618 687 !extra code for test case boundary conditions 619 688 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 620 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))689 u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 621 690 END IF 691 622 692 END DO 623 693 END DO … … 659 729 ! 660 730 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 661 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 662 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 663 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 664 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 665 & ) * 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 731 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 732 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 733 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 734 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 735 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 736 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 737 & ) / ( zbetau + 1._wp ) & 738 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 666 739 & ) * zmsk00x(ji,jj) 667 740 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 668 u_ice(ji,jj) = ( ( 669 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)670 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast671 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0672 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) 741 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 742 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 743 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 744 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 745 & ) * 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 673 746 & ) * zmsk00x(ji,jj) 674 747 ENDIF 675 748 !extra code for test case boundary conditions 676 749 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 677 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj))750 u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 678 751 END IF 679 752 END DO … … 714 787 ! 715 788 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 716 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 717 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 718 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 719 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 720 & ) * 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 789 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 790 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 791 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 792 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 793 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 794 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 795 & ) / ( zbetav + 1._wp ) & 796 & ) * 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 721 797 & ) * zmsk00y(ji,jj) 722 798 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 723 v_ice(ji,jj) = ( ( 724 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)725 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast726 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0727 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) 799 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 800 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 801 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 802 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 803 & ) * 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 728 804 & ) * zmsk00y(ji,jj) 729 805 ENDIF 730 806 !extra code for test case boundary conditions 731 807 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 732 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj))808 v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 733 809 END IF 734 810 END DO … … 744 820 ENDIF 745 821 746 !!$ IF(ln_ctl) THEN ! Convergence test 747 !!$ DO jj = 2 , jpjm1 748 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 749 !!$ END DO 750 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 751 !!$ CALL mpp_max( 'icedyn_rhg_eap', zresm ) ! max over the global domain 752 !!$ ENDIF 822 ! convergence test 823 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 753 824 ! 754 755 825 ! ! ==================== ! 756 826 END DO ! end loop over jter ! 757 827 ! ! ==================== ! 828 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 758 829 ! 759 830 CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1. ) ! only need this in ridging module after rheology completed … … 782 853 zdt2 = zdt * zdt 783 854 855 zten_i(ji,jj) = zdt 856 784 857 ! shear**2 at T points (doc eq. A16) 785 858 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 796 869 797 870 ! delta at T points 798 z delta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )799 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0800 pdelta_i(ji,jj) = z delta + rn_creepl * rswitch871 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 872 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 873 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 801 874 802 875 END DO 803 876 END DO 804 CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )805 !CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1. )877 CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 878 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. ) 806 879 807 880 ! --- Store the stress tensor for the next time step --- ! 808 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )809 881 pstress1_i (:,:) = zs1 (:,:) 810 882 pstress2_i (:,:) = zs2 (:,:) … … 815 887 ! 5) diagnostics 816 888 !------------------------------------------------------------------------------! 817 DO jj = 1, jpj818 DO ji = 1, jpi819 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice820 END DO821 END DO822 823 889 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 824 890 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 839 905 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 840 906 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 907 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 841 908 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 842 909 843 ! --- stress tensor--- !844 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN910 ! --- Stress tensor invariants (SIMIP diags) --- ! 911 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 845 912 ! 846 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )913 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 847 914 ! 848 DO jj = 2, jpjm1 849 DO ji = 2, jpim1 850 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 851 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 852 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 853 854 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 855 856 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 857 858 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 859 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 860 !! 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 861 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 862 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 863 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 864 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 915 DO jj = 1, jpj 916 DO ji = 1, jpi 917 918 ! Ice stresses 919 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 920 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 921 ! I know, this can be confusing... 922 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 923 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 924 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 925 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 926 927 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 928 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 929 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 930 865 931 END DO 866 END DO 867 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 932 END DO 868 933 ! 869 CALL iom_put( 'isig1' , zsig1 ) 870 CALL iom_put( 'isig2' , zsig2 ) 871 CALL iom_put( 'isig3' , zsig3 ) 934 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 935 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 936 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 937 938 DEALLOCATE ( zsig_I, zsig_II ) 939 940 ENDIF 941 942 ! --- Normalized stress tensor principal components --- ! 943 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 944 ! Recommendation 1 : we use ice strength, not replacement pressure 945 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 946 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 872 947 ! 873 ! Stress tensor invariants (normal and shear stress N/m) 874 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 875 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 876 877 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 948 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 949 ! 950 DO jj = 1, jpj 951 DO ji = 1, jpi 952 953 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 954 ! and **deformations** at current iterates 955 ! following Lemieux & Dupont (2020) 956 zfac = zp_delt(ji,jj) 957 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 958 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 959 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 960 961 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 962 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 963 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 964 965 ! Normalized principal stresses (used to display the ellipse) 966 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 967 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 968 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 969 END DO 970 END DO 971 ! 972 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 973 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 974 975 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 976 878 977 ENDIF 879 978 … … 949 1048 ENDIF 950 1049 ! 1050 ! --- convergence tests --- ! 1051 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 1052 IF( iom_use('uice_cvg') ) THEN 1053 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 1054 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 1055 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 1056 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 1057 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 1058 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 1059 ENDIF 1060 ENDIF 1061 ENDIF 1062 ! 1063 DEALLOCATE( zmsk00, zmsk15 ) 1064 ! 951 1065 END SUBROUTINE ice_dyn_rhg_eap 952 1066 953 SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 954 pshear, pa11, pa12, & 955 pstressp, pstressm, & 956 pstress12, strength, & 957 palphar, palphas) 1067 1068 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 1069 !!---------------------------------------------------------------------- 1070 !! *** ROUTINE rhg_cvg *** 1071 !! 1072 !! ** Purpose : check convergence of oce rheology 1073 !! 1074 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 1075 !! during the sub timestepping of rheology so as: 1076 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 1077 !! This routine is called every sub-iteration, so it is cpu expensive 1078 !! 1079 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1080 !!---------------------------------------------------------------------- 1081 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1082 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 1083 !! 1084 INTEGER :: it, idtime, istatus 1085 INTEGER :: ji, jj ! dummy loop indices 1086 REAL(wp) :: zresm ! local real 1087 CHARACTER(len=20) :: clname 1088 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 1089 !!---------------------------------------------------------------------- 1090 1091 ! create file 1092 IF( kt == nit000 .AND. kiter == 1 ) THEN 1093 ! 1094 IF( lwp ) THEN 1095 WRITE(numout,*) 1096 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 1097 WRITE(numout,*) '~~~~~~~' 1098 ENDIF 1099 ! 1100 IF( lwm ) THEN 1101 clname = 'ice_cvg.nc' 1102 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 1103 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 1104 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 1105 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 1106 istatus = NF90_ENDDEF(ncvgid) 1107 ENDIF 1108 ! 1109 ENDIF 1110 1111 ! time 1112 it = ( kt - 1 ) * kitermax + kiter 1113 1114 ! convergence 1115 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 1116 zresm = 0._wp 1117 ELSE 1118 DO jj = 1, jpj 1119 DO ji = 1, jpi 1120 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1121 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 1122 END DO 1123 END DO 1124 1125 ! cut of the boundary of the box (forced velocities) 1126 IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN 1127 zres(ji,jj) = 0._wp 1128 END IF 1129 zresm = MAXVAL( zres ) 1130 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1131 ENDIF 1132 1133 IF( lwm ) THEN 1134 ! write variables 1135 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1136 ! close file 1137 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 1138 ENDIF 1139 1140 END SUBROUTINE rhg_cvg 1141 1142 1143 SUBROUTINE update_stress_rdg( ksub, kndte, pdivu, ptension, pshear, pa11, pa12, & 1144 & pstressp, pstressm, pstress12, pstrength, palphar, palphas ) 958 1145 !!--------------------------------------------------------------------- 959 1146 !! *** SUBROUTINE update_stress_rdg *** … … 963 1150 !!--------------------------------------------------------------------- 964 1151 INTEGER, INTENT(in ) :: ksub, kndte 965 REAL(wp), INTENT(in ) :: strength1152 REAL(wp), INTENT(in ) :: pstrength 966 1153 REAL(wp), INTENT(in ) :: pdivu, ptension, pshear 967 1154 REAL(wp), INTENT(in ) :: pa11, pa12 … … 991 1178 ! Factor to maintain the same stress as in EVP (see Section 3) 992 1179 ! Can be set to 1 otherwise 993 zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 1180 ! zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 1181 zinvstressconviso = 1._wp 994 1182 995 1183 zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso … … 1009 1197 IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN 1010 1198 zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 1011 4._wp*pa12*pa12 )1199 4._wp*pa12*pa12 ) 1012 1200 1013 1201 zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma !Cos^2 … … 1084 1272 ! % range in kx 1085 1273 kx = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 1086 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1274 !!clem kx = MAX( 1, MIN( nx_yield-1, INT((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 ) ) 1275 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1087 1276 1088 1277 ky = int(zy*zinvdy) + 1 1278 !!clem ky = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) ) 1089 1279 kyw = ky - zy*zinvdy 1090 1280 1091 1281 ka = int((zatempprime-0.5_wp)*zinvda) + 1 1282 !!clem ka = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) ) 1092 1283 kaw = ka - (zatempprime-0.5_wp)*zinvda 1093 1284 1094 1285 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 1095 zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1096 & + (1._wp-zkxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1097 & + zkxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1098 & + zkxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1099 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1100 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1101 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1102 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1103 zstemp12r = zkxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1104 & + (1._wp-zkxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1105 & + zkxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1106 & + zkxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1107 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1108 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1109 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1110 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1111 zstemp22r = zkxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1112 & + (1._wp-zkxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1113 & + zkxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1114 & + zkxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1115 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1116 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1117 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1118 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1286 !!$ zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1287 !!$ & + (1._wp-zkxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1288 !!$ & + zkxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1289 !!$ & + zkxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1290 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1291 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1292 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1293 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1294 !!$ zstemp12r = zkxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1295 !!$ & + (1._wp-zkxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1296 !!$ & + zkxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1297 !!$ & + zkxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1298 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1299 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1300 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1301 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1302 !!$ zstemp22r = zkxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1303 !!$ & + (1._wp-zkxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1304 !!$ & + zkxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1305 !!$ & + zkxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1306 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1307 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1308 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1309 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1310 !!$ 1311 !!$ zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1312 !!$ & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1313 !!$ & + zkxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1314 !!$ & + zkxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1315 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1316 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1317 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1318 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1319 !!$ zstemp12s = zkxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1320 !!$ & + (1._wp-zkxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1321 !!$ & + zkxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1322 !!$ & + zkxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1323 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1324 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1325 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1326 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1327 !!$ zstemp22s = zkxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1328 !!$ & + (1._wp-zkxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1329 !!$ & + zkxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1330 !!$ & + zkxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1331 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1332 !!$ & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1333 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1334 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1335 1336 zstemp11r = s11r(kx,ky,ka) 1337 zstemp12r = s12r(kx,ky,ka) 1338 zstemp22r = s22r(kx,ky,ka) 1339 zstemp11s = s11s(kx,ky,ka) 1340 zstemp12s = s12s(kx,ky,ka) 1341 zstemp22s = s22s(kx,ky,ka) 1119 1342 1120 zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1121 & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1122 & + zkxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1123 & + zkxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1124 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1125 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1126 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1127 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1128 zstemp12s = zkxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1129 & + (1._wp-zkxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1130 & + zkxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1131 & + zkxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1132 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1133 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1134 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1135 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1136 zstemp22s = zkxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1137 & + (1._wp-zkxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1138 & + zkxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1139 & + zkxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1140 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1141 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1142 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1143 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1144 1343 1145 1344 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1146 1345 ! Tsamados 2013) 1147 1346 1148 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin1149 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin1150 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin1347 zsig11 = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin 1348 zsig12 = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin 1349 zsig22 = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1151 1350 1152 1351 !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 … … 1155 1354 1156 1355 ! Update stress 1157 zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 - 1356 zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 - 2._wp*zQ11Q12 *zsig12 1158 1357 zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12 1159 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 1358 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 2._wp*zQ11Q12 *zsig12 1160 1359 1161 1360 pstressp = zsgprm11 + zsgprm22 … … 1167 1366 IF (ksub == kndte) THEN 1168 1367 zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 1169 + zQ12Q12*zstemp22r1170 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) &1171 - zQ12Q12*zstemp12r1368 + zQ12Q12*zstemp22r 1369 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) & 1370 - zQ12Q12*zstemp12r 1172 1371 zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 1173 + zQ11Q11*zstemp22r1372 + zQ11Q11*zstemp22r 1174 1373 1175 1374 zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s & 1176 + zQ12Q12*zstemp22s1177 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) &1178 - zQ12Q12*zstemp12s1375 + zQ12Q12*zstemp22s 1376 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) & 1377 - zQ12Q12*zstemp12s 1179 1378 zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 1180 + zQ11Q11*zstemp22s1379 + zQ11Q11*zstemp22s 1181 1380 1182 1381 palphar = zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 1183 + zrotstemp22r*zdtemp221382 + zrotstemp22r*zdtemp22 1184 1383 palphas = zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 1185 + zrotstemp22s*zdtemp221384 + zrotstemp22s*zdtemp22 1186 1385 ENDIF 1187 1386 END SUBROUTINE update_stress_rdg … … 1190 1389 1191 1390 1192 SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 1193 pa11, pa12, & 1194 pmresult11, pmresult12) 1391 SUBROUTINE calc_ffrac( pstressp, pstressm, pstress12, pa11, pa12, & 1392 & pmresult11, pmresult12 ) 1195 1393 !!--------------------------------------------------------------------- 1196 1394 !! *** ROUTINE calc_ffrac *** … … 1210 1408 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1211 1409 1212 REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation 1410 !!$ REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation 1411 REAL (wp), PARAMETER :: kfrac = 1.e-3_wp ! rate of fracture formation 1213 1412 REAL (wp), PARAMETER :: threshold = 0.3_wp ! critical confinement ratio 1214 1413 !!--------------------------------------------------------------- … … 1241 1440 ENDIF 1242 1441 1243 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 & 1244 + zQ12Q12*zsigma22 ! S(1,1) 1245 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 & 1246 + zQ11Q11*zsigma22 ! S(2,2) 1442 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 + zQ12Q12*zsigma22 ! S(1,1) 1443 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 + zQ11Q11*zsigma22 ! S(2,2) 1247 1444 1248 1445 ! Pure divergence … … 1255 1452 ! which leads to the loss of their shape, so we again model it through diffusion 1256 1453 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1257 pmresult11 = kfrac * (pa11 - zQ12Q12)1258 pmresult12 = kfrac * (pa12 + zQ11Q12)1454 pmresult11 = - kfrac * (pa11 - zQ12Q12) 1455 pmresult12 = - kfrac * (pa12 + zQ11Q12) 1259 1456 1260 1457 ! Shear faulting … … 1263 1460 pmresult12 = 0.0_wp 1264 1461 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 1265 pmresult11 = kfrac * (pa11 - zQ12Q12)1266 pmresult12 = kfrac * (pa12 + zQ11Q12)1462 pmresult11 = - kfrac * (pa11 - zQ12Q12) 1463 pmresult12 = - kfrac * (pa12 + zQ11Q12) 1267 1464 1268 1465 ! Horizontal spalling … … 1295 1492 REAL(wp) :: da, dx, dy, dp, dz, a1 1296 1493 1494 !!clem 1495 REAL(wp) :: zw1, zw2, zfac, ztemp 1496 REAL(wp) :: idx, idy, idz 1497 1297 1498 REAL(wp), PARAMETER :: eps6 = 1.0e-6_wp 1298 1499 !!---------------------------------------------------------------------- … … 1305 1506 id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. ) 1306 1507 id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. ) 1307 id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. )1308 id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. )1508 id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. ) 1509 id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. ) 1309 1510 ! 1310 1511 IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN ! fields exist … … 1312 1513 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i ) 1313 1514 CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 1314 CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11)1315 CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12)1515 CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11 ) 1516 CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12 ) 1316 1517 ELSE ! start rheology from rest 1317 1518 IF(lwp) WRITE(numout,*) … … 1320 1521 stress2_i (:,:) = 0._wp 1321 1522 stress12_i(:,:) = 0._wp 1322 aniso_11 (:,:) = 0.5_wp1323 aniso_12 (:,:) = 0._wp1523 aniso_11 (:,:) = 0.5_wp 1524 aniso_12 (:,:) = 0._wp 1324 1525 ENDIF 1325 1526 ELSE !* Start from rest … … 1329 1530 stress2_i (:,:) = 0._wp 1330 1531 stress12_i(:,:) = 0._wp 1331 aniso_11 (:,:) = 0.5_wp1332 aniso_12 (:,:) = 0._wp1532 aniso_11 (:,:) = 0.5_wp 1533 aniso_12 (:,:) = 0._wp 1333 1534 ENDIF 1334 1535 ! 1335 1536 1336 da = 0.5_wp/real(na_yield-1,kind=wp) 1337 ainit = 0.5_wp - da 1338 dx = rpi/real(nx_yield-1,kind=wp) 1339 xinit = rpi + 0.25_wp*rpi - dx 1340 dz = rpi/real(nz,kind=wp) 1341 zinit = -rpi*0.5_wp 1342 dy = rpi/real(ny_yield-1,kind=wp) 1343 yinit = -dy 1344 1345 DO ia=1,na_yield 1346 DO ix=1,nx_yield 1347 DO iy=1,ny_yield 1348 s11r(ix,iy,ia) = 0._wp 1349 s12r(ix,iy,ia) = 0._wp 1350 s22r(ix,iy,ia) = 0._wp 1351 s11s(ix,iy,ia) = 0._wp 1352 s12s(ix,iy,ia) = 0._wp 1353 s22s(ix,iy,ia) = 0._wp 1354 IF (ia <= na_yield-1) THEN 1355 DO iz=1,nz 1356 s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1357 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1358 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1359 s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1360 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1361 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1362 s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1363 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1364 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1365 s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1366 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1367 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1368 s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1369 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1370 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1371 s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1372 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1373 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1374 ENDDO 1375 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1376 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1377 IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1378 IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1379 IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1380 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1381 ELSE 1382 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1383 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1384 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1385 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1386 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1387 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1388 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1389 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1390 IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1391 IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1392 IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1393 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1394 ENDIF 1395 ENDDO 1396 ENDDO 1397 ENDDO 1537 da = 0.5_wp/real(na_yield-1,kind=wp) 1538 ainit = 0.5_wp - da 1539 dx = rpi/real(nx_yield-1,kind=wp) 1540 xinit = rpi + 0.25_wp*rpi - dx 1541 dz = rpi/real(nz,kind=wp) 1542 zinit = -rpi*0.5_wp 1543 dy = rpi/real(ny_yield-1,kind=wp) 1544 yinit = -dy 1545 1546 s11r(:,:,:) = 0._wp 1547 s12r(:,:,:) = 0._wp 1548 s22r(:,:,:) = 0._wp 1549 s11s(:,:,:) = 0._wp 1550 s12s(:,:,:) = 0._wp 1551 s22s(:,:,:) = 0._wp 1552 1553 !!$ DO ia=1,na_yield 1554 !!$ DO ix=1,nx_yield 1555 !!$ DO iy=1,ny_yield 1556 !!$ s11r(ix,iy,ia) = 0._wp 1557 !!$ s12r(ix,iy,ia) = 0._wp 1558 !!$ s22r(ix,iy,ia) = 0._wp 1559 !!$ s11s(ix,iy,ia) = 0._wp 1560 !!$ s12s(ix,iy,ia) = 0._wp 1561 !!$ s22s(ix,iy,ia) = 0._wp 1562 !!$ IF (ia <= na_yield-1) THEN 1563 !!$ DO iz=1,nz 1564 !!$ s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1565 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1566 !!$ s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1567 !!$ s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1568 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1569 !!$ s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1570 !!$ s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1571 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1572 !!$ s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1573 !!$ s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1574 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1575 !!$ s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1576 !!$ s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1577 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1578 !!$ s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1579 !!$ s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1580 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1581 !!$ s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1582 !!$ ENDDO 1583 !!$ IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1584 !!$ IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1585 !!$ IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1586 !!$ IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1587 !!$ IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1588 !!$ IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1589 !!$ ELSE 1590 !!$ s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1591 !!$ s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1592 !!$ s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1593 !!$ s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1594 !!$ s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1595 !!$ s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1596 !!$ IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1597 !!$ IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 1598 !!$ IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 1599 !!$ IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 1600 !!$ IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 1601 !!$ IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1602 !!$ ENDIF 1603 !!$ ENDDO 1604 !!$ ENDDO 1605 !!$ ENDDO 1606 1607 !! faster but still very slow => to be improved 1608 zfac = dz/sin(2._wp*pphi) 1609 DO ia = 1, na_yield-1 1610 zw1 = w1(ainit+ia*da) 1611 zw2 = w2(ainit+ia*da) 1612 DO iz = 1, nz 1613 idz = zinit+iz*dz 1614 ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 1615 DO iy = 1, ny_yield 1616 idy = yinit+iy*dy 1617 DO ix = 1, nx_yield 1618 idx = xinit+ix*dx 1619 s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac 1620 s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac 1621 s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac 1622 s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac 1623 s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac 1624 s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac 1625 END DO 1626 END DO 1627 END DO 1628 END DO 1629 1630 zfac = 1._wp/sin(2._wp*pphi) 1631 ia = na_yield 1632 DO iy = 1, ny_yield 1633 idy = yinit+iy*dy 1634 DO ix = 1, nx_yield 1635 idx = xinit+ix*dx 1636 s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac 1637 s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac 1638 s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac 1639 s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac 1640 s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac 1641 s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac 1642 ENDDO 1643 ENDDO 1644 WHERE (ABS(s11r(:,:,:)) < eps6) s11r(:,:,:) = 0._wp 1645 WHERE (ABS(s12r(:,:,:)) < eps6) s12r(:,:,:) = 0._wp 1646 WHERE (ABS(s22r(:,:,:)) < eps6) s22r(:,:,:) = 0._wp 1647 WHERE (ABS(s11s(:,:,:)) < eps6) s11s(:,:,:) = 0._wp 1648 WHERE (ABS(s12s(:,:,:)) < eps6) s12s(:,:,:) = 0._wp 1649 WHERE (ABS(s22s(:,:,:)) < eps6) s22s(:,:,:) = 0._wp 1650 1398 1651 1399 1652 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file … … 1405 1658 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) 1406 1659 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) 1407 CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11)1408 CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12)1660 CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11 ) 1661 CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12 ) 1409 1662 ! 1410 1663 ENDIF … … 1421 1674 !!------------------------------------------------------------------- 1422 1675 1423 w1 = - 223.87569446_wp &1424 & + 2361.2198663_wp*pa &1676 w1 = - 223.87569446_wp & 1677 & + 2361.21986630_wp*pa & 1425 1678 & - 10606.56079975_wp*pa*pa & 1426 1679 & + 26315.50025642_wp*pa*pa*pa & … … 1428 1681 & + 34397.72407466_wp*pa*pa*pa*pa*pa & 1429 1682 & - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 1430 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa1683 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 1431 1684 1432 1685 END FUNCTION w1 … … 1441 1694 !!------------------------------------------------------------------- 1442 1695 1443 w2 = - 6670.68911883_wp &1444 & + 70222.33061536_wp*pa &1445 & - 314871.71525448_wp*pa*pa &1446 & + 779570.02793492_wp*pa*pa*pa &1696 w2 = - 6670.68911883_wp & 1697 & + 70222.33061536_wp*pa & 1698 & - 314871.71525448_wp*pa*pa & 1699 & + 779570.02793492_wp*pa*pa*pa & 1447 1700 & - 1151098.82436864_wp*pa*pa*pa*pa & 1448 1701 & + 1013896.59464498_wp*pa*pa*pa*pa*pa & 1449 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa &1450 & + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa1702 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 1703 & + 102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa 1451 1704 1452 1705 END FUNCTION w2 … … 1806 2059 !! Default option Empty module NO SI3 sea-ice model 1807 2060 !!---------------------------------------------------------------------- 1808 1809 2061 #endif 2062 1810 2063 !!============================================================================== 1811 2064 END MODULE icedyn_rhg_eap
Note: See TracChangeset
for help on using the changeset viewer.
