Changeset 13472 for NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2020-09-16T15:05:19+02:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r13461 r13472 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 50 51 # include "do_loop_substitute.h90" 51 52 # include "domzgr_substitute.h90" 53 54 !! for convergence tests 55 INTEGER :: ncvgid ! netcdf file id 56 INTEGER :: nvarid ! netcdf variable id 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 52 58 !!---------------------------------------------------------------------- 53 59 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 121 127 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 122 128 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 129 REAl(wp) :: zbetau, zbetav 123 130 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 124 131 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars … … 127 134 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 128 135 ! 129 REAL(wp) :: zresm ! Maximal error on ice velocity130 136 REAL(wp) :: zintb, zintn ! dummy argument 131 137 REAL(wp) :: zfac_x, zfac_y … … 143 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 144 150 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 145 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence146 151 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 147 152 ! ! ocean surface (ssh_m) if ice is not embedded … … 162 167 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 163 168 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 169 !! --- check convergence 170 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 164 171 !! --- diags 165 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00166 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 167 173 !! --- SIMIP diags … … 176 182 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 177 183 ! 178 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 184 ! for diagnostics and convergence tests 185 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 186 DO_2D( 1, 1, 1, 1 ) 187 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 188 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 189 END_2D 190 ! 191 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 179 192 !------------------------------------------------------------------------------! 180 193 ! 0) mask at F points for the ice … … 220 233 z1_ecc2 = 1._wp / ecc2 221 234 222 ! Time step for subcycling223 zdtevp = rDt_ice / REAL( nn_nevp )224 z1_dtevp = 1._wp / zdtevp225 226 235 ! alpha parameters (Bouillon 2009) 227 236 IF( .NOT. ln_aEVP ) THEN 228 zalph1 = ( 2._wp * rn_relast * rDt_ice ) * z1_dtevp 237 zdtevp = rDt_ice / REAL( nn_nevp ) 238 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 229 239 zalph2 = zalph1 * z1_ecc2 230 240 231 241 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 232 242 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 243 ELSE 244 zdtevp = rdt_ice 245 ! zalpha parameters set later on adaptatively 233 246 ENDIF 247 z1_dtevp = 1._wp / zdtevp 234 248 235 249 ! Initialise stress tensor … … 242 256 243 257 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 244 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile258 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 245 259 ELSE ; zkt = 0._wp 246 260 ENDIF … … 310 324 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 311 325 ! ice-bottom stress at U points 312 zvCr = zaU(ji,jj) * rn_ depfra * hu(ji,jj,Kmm)313 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )326 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 327 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 314 328 ! ice-bottom stress at V points 315 zvCr = zaV(ji,jj) * rn_ depfra * hv(ji,jj,Kmm)316 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )329 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 330 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 317 331 ! ice_bottom stress at T points 318 zvCr = at_i(ji,jj) * rn_ depfra * ht(ji,jj)319 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )332 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 333 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 320 334 END_2D 321 335 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp ) … … 337 351 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 338 352 ! 339 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 340 !!$ DO jj = 1, jpjm1 341 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 342 !!$ zv_ice(:,jj) = v_ice(:,jj) 343 !!$ END DO 344 !!$ ENDIF 353 ! convergence test 354 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 355 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 356 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 357 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 358 END_2D 359 ENDIF 345 360 346 361 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 380 395 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 381 396 382 ! alpha & betafor aEVP397 ! alpha for aEVP 383 398 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 384 399 ! alpha = beta = sqrt(4*gamma) … … 388 403 zalph2 = zalph1 389 404 z1_alph2 = z1_alph1 405 ! explicit: 406 ! z1_alph1 = 1._wp / zalph1 407 ! z1_alph2 = 1._wp / zalph1 408 ! zalph1 = zalph1 - 1._wp 409 ! zalph2 = zalph1 390 410 ENDIF 391 411 … … 397 417 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 398 418 419 ! Save beta at T-points for further computations 420 IF( ln_aEVP ) THEN 421 DO_2D( 1, 1, 1, 1 ) 422 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 423 END_2D 424 ENDIF 425 399 426 DO_2D( 1, 0, 1, 0 ) 400 427 401 ! alpha & betafor aEVP428 ! alpha for aEVP 402 429 IF( ln_aEVP ) THEN 403 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj)) )430 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 404 431 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 405 zbeta(ji,jj) = zalph2 432 ! explicit: 433 ! z1_alph2 = 1._wp / zalph2 434 ! zalph2 = zalph2 - 1._wp 406 435 ENDIF 407 436 … … 469 498 ! 470 499 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 471 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 472 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 473 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 474 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 475 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 500 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 501 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 502 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 503 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 504 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 505 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 506 & ) / ( zbetav + 1._wp ) & 507 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 476 508 & ) * zmsk00y(ji,jj) 477 509 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 478 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity479 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)480 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast481 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0482 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin483 & ) 510 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 511 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 512 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 513 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 514 & ) * 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 515 & ) * zmsk00y(ji,jj) 484 516 ENDIF 485 517 END_2D … … 518 550 ! 519 551 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 520 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 521 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 522 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 523 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 524 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 552 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 553 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 554 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 555 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 556 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 557 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 558 & ) / ( zbetau + 1._wp ) & 559 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 525 560 & ) * zmsk00x(ji,jj) 526 561 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 527 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity528 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)529 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast530 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0531 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin532 & 562 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 563 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 564 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 565 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 566 & ) * 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 567 & ) * zmsk00x(ji,jj) 533 568 ENDIF 534 569 END_2D … … 569 604 ! 570 605 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 571 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 572 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 573 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 574 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 575 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 606 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 607 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 608 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 609 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 610 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 611 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 612 & ) / ( zbetau + 1._wp ) & 613 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 576 614 & ) * zmsk00x(ji,jj) 577 615 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 578 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity579 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)580 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast581 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0582 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin583 & 616 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 617 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 618 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 619 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 620 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 621 & ) * zmsk00x(ji,jj) 584 622 ENDIF 585 623 END_2D … … 618 656 ! 619 657 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 620 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 621 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 622 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 623 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 624 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 658 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 659 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 660 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 661 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 662 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 663 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 664 & ) / ( zbetav + 1._wp ) & 665 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 625 666 & ) * zmsk00y(ji,jj) 626 667 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 627 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity628 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)629 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast630 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0631 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin632 & 668 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 669 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 670 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 671 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 672 & ) * 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 673 & ) * zmsk00y(ji,jj) 633 674 ENDIF 634 675 END_2D … … 643 684 ENDIF 644 685 645 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 646 !!$ DO jj = 2 , jpjm1 647 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 648 !!$ END DO 649 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 650 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 651 !!$ ENDIF 686 ! convergence test 687 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 652 688 ! 653 689 ! ! ==================== ! 654 690 END DO ! end loop over jter ! 655 691 ! ! ==================== ! 692 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 656 693 ! 657 694 !------------------------------------------------------------------------------! … … 706 743 ! 5) diagnostics 707 744 !------------------------------------------------------------------------------! 708 DO_2D( 1, 1, 1, 1 )709 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice710 END_2D711 712 745 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 713 746 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 764 797 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 765 798 ENDIF 766 799 767 800 ! --- SIMIP --- ! 768 801 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 818 851 ENDIF 819 852 ! 853 ! --- convergence tests --- ! 854 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 855 IF( iom_use('uice_cvg') ) THEN 856 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 857 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 858 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 859 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 860 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 861 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 862 ENDIF 863 ENDIF 864 ENDIF 865 ! 866 DEALLOCATE( zmsk00, zmsk15 ) 867 ! 820 868 END SUBROUTINE ice_dyn_rhg_evp 869 870 871 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 872 !!---------------------------------------------------------------------- 873 !! *** ROUTINE rhg_cvg *** 874 !! 875 !! ** Purpose : check convergence of oce rheology 876 !! 877 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 878 !! during the sub timestepping of rheology so as: 879 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 880 !! This routine is called every sub-iteration, so it is cpu expensive 881 !! 882 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 883 !!---------------------------------------------------------------------- 884 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 885 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 886 !! 887 INTEGER :: it, idtime, istatus 888 INTEGER :: ji, jj ! dummy loop indices 889 REAL(wp) :: zresm ! local real 890 CHARACTER(len=20) :: clname 891 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 892 !!---------------------------------------------------------------------- 893 894 ! create file 895 IF( kt == nit000 .AND. kiter == 1 ) THEN 896 ! 897 IF( lwp ) THEN 898 WRITE(numout,*) 899 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 900 WRITE(numout,*) '~~~~~~~' 901 ENDIF 902 ! 903 IF( lwm ) THEN 904 clname = 'ice_cvg.nc' 905 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 906 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 907 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 908 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 909 istatus = NF90_ENDDEF(ncvgid) 910 ENDIF 911 ! 912 ENDIF 913 914 ! time 915 it = ( kt - 1 ) * kitermax + kiter 916 917 ! convergence 918 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 919 zresm = 0._wp 920 ELSE 921 DO_2D( 1, 1, 1, 1 ) 922 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 923 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 924 END_2D 925 zresm = MAXVAL( zres ) 926 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 927 ENDIF 928 929 IF( lwm ) THEN 930 ! write variables 931 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 932 ! close file 933 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 934 ENDIF 935 936 END SUBROUTINE rhg_cvg 821 937 822 938 … … 876 992 END SUBROUTINE rhg_evp_rst 877 993 994 878 995 #else 879 996 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
