- Timestamp:
- 2020-09-15T09:27:47+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90
r13271 r13466 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 49 50 !! * Substitutions 50 51 # include "vectopt_loop_substitute.h90" 52 53 !! for convergence tests 54 INTEGER :: ncvgid ! netcdf file id 55 INTEGER :: nvarid ! netcdf variable id 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 51 57 !!---------------------------------------------------------------------- 52 58 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 119 125 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 120 126 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 127 REAl(wp) :: zbetau, zbetav 121 128 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 122 129 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars … … 125 132 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 126 133 ! 127 REAL(wp) :: zresm ! Maximal error on ice velocity128 134 REAL(wp) :: zintb, zintn ! dummy argument 129 135 REAL(wp) :: zfac_x, zfac_y … … 141 147 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 142 148 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 143 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence144 149 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 145 150 ! ! ocean surface (ssh_m) if ice is not embedded … … 160 165 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 161 166 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 167 !! --- check convergence 168 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 162 169 !! --- diags 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00164 170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 165 171 !! --- SIMIP diags … … 174 180 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 175 181 ! 176 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 182 ! for diagnostics and convergence tests 183 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 184 DO jj = 1, jpj 185 DO ji = 1, jpi 186 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 187 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 188 END DO 189 END DO 190 ! 191 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 177 192 !------------------------------------------------------------------------------! 178 193 ! 0) mask at F points for the ice … … 222 237 z1_ecc2 = 1._wp / ecc2 223 238 224 ! Time step for subcycling225 zdtevp = rdt_ice / REAL( nn_nevp )226 z1_dtevp = 1._wp / zdtevp227 228 239 ! alpha parameters (Bouillon 2009) 229 240 IF( .NOT. ln_aEVP ) THEN 230 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 241 zdtevp = rdt_ice / REAL( nn_nevp ) 242 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 231 243 zalph2 = zalph1 * z1_ecc2 232 244 233 245 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 234 246 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 247 ELSE 248 zdtevp = rdt_ice 249 ! zalpha parameters set later on adaptatively 235 250 ENDIF 251 z1_dtevp = 1._wp / zdtevp 236 252 237 253 ! Initialise stress tensor … … 244 260 245 261 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 246 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile262 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 247 263 ELSE ; zkt = 0._wp 248 264 ENDIF … … 315 331 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) 316 332 ! ice-bottom stress at U points 317 zvCr = zaU(ji,jj) * rn_ depfra * hu_n(ji,jj)318 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )333 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 334 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 319 335 ! ice-bottom stress at V points 320 zvCr = zaV(ji,jj) * rn_ depfra * hv_n(ji,jj)321 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )336 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 337 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 322 338 ! ice_bottom stress at T points 323 zvCr = at_i(ji,jj) * rn_ depfra * ht_n(ji,jj)324 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )339 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 340 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 325 341 END DO 326 342 END DO … … 345 361 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 346 362 ! 347 !!$ IF(ln_ctl) THEN ! Convergence test 348 !!$ DO jj = 1, jpjm1 349 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 350 !!$ zv_ice(:,jj) = v_ice(:,jj) 351 !!$ END DO 352 !!$ ENDIF 363 ! convergence test 364 IF( ln_rhg_chkcvg ) THEN 365 DO jj = 1, jpj 366 DO ji = 1, jpi 367 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 368 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 369 END DO 370 END DO 371 ENDIF 353 372 354 373 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 391 410 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 392 411 393 ! alpha & betafor aEVP412 ! alpha for aEVP 394 413 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 395 414 ! alpha = beta = sqrt(4*gamma) … … 399 418 zalph2 = zalph1 400 419 z1_alph2 = z1_alph1 420 ! explicit: 421 ! z1_alph1 = 1._wp / zalph1 422 ! z1_alph2 = 1._wp / zalph1 423 ! zalph1 = zalph1 - 1._wp 424 ! zalph2 = zalph1 401 425 ENDIF 402 426 … … 409 433 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 410 434 435 ! Save beta at T-points for further computations 436 IF( ln_aEVP ) THEN 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 440 END DO 441 END DO 442 ENDIF 443 411 444 DO jj = 1, jpjm1 412 445 DO ji = 1, jpim1 413 446 414 ! alpha & betafor aEVP447 ! alpha for aEVP 415 448 IF( ln_aEVP ) THEN 416 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj)) )449 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 417 450 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 418 zbeta(ji,jj) = zalph2 451 ! explicit: 452 ! z1_alph2 = 1._wp / zalph2 453 ! zalph2 = zalph2 - 1._wp 419 454 ENDIF 420 455 … … 486 521 ! 487 522 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 488 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 489 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 490 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 491 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 492 & ) * 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 523 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 524 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 525 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 526 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 527 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 528 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 529 & ) / ( zbetav + 1._wp ) & 530 & ) * 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 493 531 & ) * zmsk00y(ji,jj) 494 532 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 495 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity496 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)497 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast498 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0499 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin500 & ) 533 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 534 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 535 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 536 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 537 & ) * 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 538 & ) * zmsk00y(ji,jj) 501 539 ENDIF 502 540 END DO … … 537 575 ! 538 576 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 539 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 540 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 541 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 542 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 543 & ) * 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 577 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 578 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 579 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 580 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 581 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 582 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 583 & ) / ( zbetau + 1._wp ) & 584 & ) * 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 544 585 & ) * zmsk00x(ji,jj) 545 586 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 546 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity547 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)548 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast549 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0550 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin551 & 587 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 588 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 589 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 590 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 591 & ) * 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 592 & ) * zmsk00x(ji,jj) 552 593 ENDIF 553 594 END DO … … 590 631 ! 591 632 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 592 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 593 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 594 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 595 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 596 & ) * 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 633 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 634 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 635 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 636 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 637 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 638 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 639 & ) / ( zbetau + 1._wp ) & 640 & ) * 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 597 641 & ) * zmsk00x(ji,jj) 598 642 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 599 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity600 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)601 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast602 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0603 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin604 & 643 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 644 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 645 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 646 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 647 & ) * 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 648 & ) * zmsk00x(ji,jj) 605 649 ENDIF 606 650 END DO … … 641 685 ! 642 686 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 643 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 644 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 645 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 646 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 647 & ) * 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 687 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 688 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 689 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 690 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 691 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 692 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 693 & ) / ( zbetav + 1._wp ) & 694 & ) * 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 648 695 & ) * zmsk00y(ji,jj) 649 696 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 650 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity651 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)652 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast653 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0654 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin655 & 697 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 698 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 699 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 700 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 701 & ) * 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 702 & ) * zmsk00y(ji,jj) 656 703 ENDIF 657 704 END DO … … 667 714 ENDIF 668 715 669 !!$ IF(ln_ctl) THEN ! Convergence test 670 !!$ DO jj = 2 , jpjm1 671 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 672 !!$ END DO 673 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 674 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 675 !!$ ENDIF 716 ! convergence test 717 IF( ln_rhg_chkcvg ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 676 718 ! 677 719 ! ! ==================== ! 678 720 END DO ! end loop over jter ! 679 721 ! ! ==================== ! 722 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 680 723 ! 681 724 !------------------------------------------------------------------------------! … … 734 777 ! 5) diagnostics 735 778 !------------------------------------------------------------------------------! 736 DO jj = 1, jpj737 DO ji = 1, jpi738 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice739 END DO740 END DO741 742 779 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 743 780 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 796 833 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 797 834 ENDIF 798 835 799 836 ! --- SIMIP --- ! 800 837 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 852 889 ENDIF 853 890 ! 891 ! --- convergence tests --- ! 892 IF( ln_rhg_chkcvg ) THEN 893 IF( iom_use('uice_cvg') ) THEN 894 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 895 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 896 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 897 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 898 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 899 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 900 ENDIF 901 ENDIF 902 ENDIF 903 ! 904 DEALLOCATE( zmsk00, zmsk15 ) 905 ! 854 906 END SUBROUTINE ice_dyn_rhg_evp 907 908 909 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 910 !!---------------------------------------------------------------------- 911 !! *** ROUTINE rhg_cvg *** 912 !! 913 !! ** Purpose : check convergence of oce rheology 914 !! 915 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 916 !! during the sub timestepping of rheology so as: 917 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 918 !! This routine is called every sub-iteration, so it is cpu expensive 919 !! 920 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 921 !!---------------------------------------------------------------------- 922 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 923 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 924 !! 925 INTEGER :: it, idtime, istatus 926 INTEGER :: ji, jj ! dummy loop indices 927 REAL(wp) :: zresm ! local real 928 CHARACTER(len=20) :: clname 929 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 930 !!---------------------------------------------------------------------- 931 932 ! create file 933 IF( kt == nit000 .AND. kiter == 1 ) THEN 934 ! 935 IF( lwp ) THEN 936 WRITE(numout,*) 937 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 938 WRITE(numout,*) '~~~~~~~' 939 ENDIF 940 ! 941 IF( lwm ) THEN 942 clname = 'ice_cvg.nc' 943 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 944 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 945 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 946 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 947 istatus = NF90_ENDDEF(ncvgid) 948 ENDIF 949 ! 950 ENDIF 951 952 ! time 953 it = ( kt - 1 ) * kitermax + kiter 954 955 ! convergence 956 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 957 zresm = 0._wp 958 ELSE 959 DO jj = 1, jpj 960 DO ji = 1, jpi 961 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 962 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 963 END DO 964 END DO 965 zresm = MAXVAL( zres ) 966 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 967 ENDIF 968 969 IF( lwm ) THEN 970 ! write variables 971 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 972 ! close file 973 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 974 ENDIF 975 976 END SUBROUTINE rhg_cvg 855 977 856 978 … … 910 1032 END SUBROUTINE rhg_evp_rst 911 1033 1034 912 1035 #else 913 1036 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.