- Timestamp:
- 2020-10-06T18:17:44+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90
r13570 r13571 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 50 51 # include "do_loop_substitute.h90" 51 52 # include "domzgr_substitute.h90" 53 54 !! for convergence tests 55 INTEGER :: ncvgid ! netcdf file id 56 INTEGER :: nvarid ! netcdf variable id 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 52 58 !!---------------------------------------------------------------------- 53 59 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 121 127 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 122 128 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 129 REAl(wp) :: zbetau, zbetav 123 130 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 124 REAL(wp) :: z delta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars131 REAL(wp) :: zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 125 132 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 126 133 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 127 134 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 128 135 ! 129 REAL(wp) :: zresm ! Maximal error on ice velocity130 136 REAL(wp) :: zintb, zintn ! dummy argument 131 137 REAL(wp) :: zfac_x, zfac_y 132 138 REAL(wp) :: zshear, zdum1, zdum2 133 139 ! 134 REAL(wp), DIMENSION(jpi,jpj) :: z p_delt !P/delta at T points140 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 135 141 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 136 142 ! … … 139 145 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 140 146 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 141 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 147 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 142 148 ! 143 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 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 … … 157 162 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 158 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 159 REAL(wp), DIMENSION(jpi,jpj) :: zfmask , zwf! mask at F points for the ice164 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 160 165 161 166 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 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 … … 191 204 192 205 ! Lateral boundary conditions on velocity (modify zfmask) 193 zwf(:,:) = zfmask(:,:)194 206 DO_2D( 0, 0, 0, 0 ) 195 207 IF( zfmask(ji,jj) == 0._wp ) THEN 196 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) ) ) 208 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 209 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 197 210 ENDIF 198 211 END_2D 199 212 DO jj = 2, jpjm1 200 213 IF( zfmask(1,jj) == 0._wp ) THEN 201 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )214 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 202 215 ENDIF 203 216 IF( zfmask(jpi,jj) == 0._wp ) THEN 204 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )205 217 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 218 ENDIF 206 219 END DO 207 220 DO ji = 2, jpim1 208 221 IF( zfmask(ji,1) == 0._wp ) THEN 209 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )222 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 210 223 ENDIF 211 224 IF( zfmask(ji,jpj) == 0._wp ) THEN 212 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )225 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 213 226 ENDIF 214 227 END DO … … 228 241 z1_ecc2 = 1._wp / ecc2 229 242 230 ! Time step for subcycling231 zdtevp = rDt_ice / REAL( nn_nevp )232 z1_dtevp = 1._wp / zdtevp233 234 243 ! alpha parameters (Bouillon 2009) 235 244 IF( .NOT. ln_aEVP ) THEN 236 zalph1 = ( 2._wp * rn_relast * rDt_ice ) * z1_dtevp 245 zdtevp = rDt_ice / REAL( nn_nevp ) 246 zalph1 = 2._wp * rn_relast * REAL( nn_nevp ) 237 247 zalph2 = zalph1 * z1_ecc2 238 248 239 249 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 240 250 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 251 ELSE 252 zdtevp = rdt_ice 253 ! zalpha parameters set later on adaptatively 241 254 ENDIF 255 z1_dtevp = 1._wp / zdtevp 242 256 243 257 ! Initialise stress tensor … … 250 264 251 265 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 252 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile266 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 253 267 ELSE ; zkt = 0._wp 254 268 ENDIF … … 322 336 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) 323 337 ! ice-bottom stress at U points 324 zvCr = zaU(ji,jj) * rn_ depfra * hu(ji,jj,Kmm)325 ztaux_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )338 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 339 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 326 340 ! ice-bottom stress at V points 327 zvCr = zaV(ji,jj) * rn_ depfra * hv(ji,jj,Kmm)328 ztauy_base(ji,jj) = - rn_ icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )341 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 342 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 329 343 ! ice_bottom stress at T points 330 zvCr = at_i(ji,jj) * rn_ depfra * ht(ji,jj)331 tau_icebfr(ji,jj) = - rn_ icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )344 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 345 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 332 346 END_2D 333 347 #if defined key_mpi3 … … 353 367 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 354 368 ! 355 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 356 !!$ DO jj = 1, jpjm1 357 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 358 !!$ zv_ice(:,jj) = v_ice(:,jj) 359 !!$ END DO 360 !!$ ENDIF 369 ! convergence test 370 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 371 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 372 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 373 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 374 END_2D 375 ENDIF 361 376 362 377 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 369 384 370 385 END_2D 371 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 372 373 DO_2D( 0, 1, 0, 1 ) 386 387 DO_2D( 0, 0, 0, 0 ) 374 388 375 389 ! shear**2 at T points (doc eq. A16) … … 391 405 392 406 ! delta at T points 393 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 394 395 ! P/delta at T points 396 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 397 398 ! alpha & beta for aEVP 407 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 408 409 END_2D 410 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 411 412 ! P/delta at T points 413 DO_2D( 1, 1, 1, 1 ) 414 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 415 END_2D 416 417 DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 418 419 ! divergence at T points (duplication to avoid communications) 420 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 421 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 422 & ) * r1_e1e2t(ji,jj) 423 424 ! tension at T points (duplication to avoid communications) 425 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) & 426 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 427 & ) * r1_e1e2t(ji,jj) 428 429 ! alpha for aEVP 399 430 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 400 431 ! alpha = beta = sqrt(4*gamma) … … 404 435 zalph2 = zalph1 405 436 z1_alph2 = z1_alph1 437 ! explicit: 438 ! z1_alph1 = 1._wp / zalph1 439 ! z1_alph2 = 1._wp / zalph1 440 ! zalph1 = zalph1 - 1._wp 441 ! zalph2 = zalph1 406 442 ENDIF 407 443 408 444 ! stress at T points (zkt/=0 if landfast) 409 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta *(1._wp - zkt) ) ) * z1_alph1410 zs2(ji,jj) = ( zs2(ji,jj) *zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2445 zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 446 zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 411 447 412 448 END_2D 413 #if defined key_mpi3 414 CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 415 #else 416 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 417 #endif 449 ! Save beta at T-points for further computations 450 IF( ln_aEVP ) THEN 451 DO_2D( 1, 1, 1, 1 ) 452 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 453 END_2D 454 ENDIF 455 418 456 DO_2D( 1, 0, 1, 0 ) 419 457 420 ! alpha & betafor aEVP458 ! alpha for aEVP 421 459 IF( ln_aEVP ) THEN 422 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj)) )460 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 423 461 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 424 zbeta(ji,jj) = zalph2 462 ! explicit: 463 ! z1_alph2 = 1._wp / zalph2 464 ! zalph2 = zalph2 - 1._wp 425 465 ENDIF 426 466 … … 488 528 ! 489 529 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 490 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 491 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 492 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 493 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 494 & ) * 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 530 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 531 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 532 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 533 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 534 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 535 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 536 & ) / ( zbetav + 1._wp ) & 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 495 538 & ) * zmsk00y(ji,jj) 496 539 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 497 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity498 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)499 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast500 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0501 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin502 & ) 540 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 541 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 542 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 543 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 544 & ) * 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 545 & ) * zmsk00y(ji,jj) 503 546 ENDIF 504 547 END_2D … … 541 584 ! 542 585 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 543 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 544 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 545 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 546 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 547 & ) * 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 586 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 587 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(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) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 590 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 591 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 592 & ) / ( zbetau + 1._wp ) & 593 & ) * 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 548 594 & ) * zmsk00x(ji,jj) 549 595 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 550 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity551 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)552 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast553 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0554 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin555 & 596 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 597 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 598 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 599 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 600 & ) * 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 601 & ) * zmsk00x(ji,jj) 556 602 ENDIF 557 603 END_2D … … 596 642 ! 597 643 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 598 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 599 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 600 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 601 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 602 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 644 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 645 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 646 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 647 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 648 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 649 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 650 & ) / ( zbetau + 1._wp ) & 651 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 603 652 & ) * zmsk00x(ji,jj) 604 653 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 605 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity606 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)607 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast608 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0609 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin610 & 654 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 655 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 656 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 657 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 658 & ) * 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 659 & ) * zmsk00x(ji,jj) 611 660 ENDIF 612 661 END_2D … … 649 698 ! 650 699 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 651 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 652 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 653 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 654 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 655 & ) * 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 700 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 701 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 702 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 703 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 704 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 705 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 706 & ) / ( zbetav + 1._wp ) & 707 & ) * 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 656 708 & ) * zmsk00y(ji,jj) 657 709 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 658 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity659 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)660 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast661 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0662 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin663 & 710 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 711 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 712 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 713 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 714 & ) * 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 715 & ) * zmsk00y(ji,jj) 664 716 ENDIF 665 717 END_2D … … 678 730 ENDIF 679 731 680 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 681 !!$ DO jj = 2 , jpjm1 682 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 683 !!$ END DO 684 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 685 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 686 !!$ ENDIF 732 ! convergence test 733 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 687 734 ! 688 735 ! ! ==================== ! 689 736 END DO ! end loop over jter ! 690 737 ! ! ==================== ! 738 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 691 739 ! 692 740 !------------------------------------------------------------------------------! … … 702 750 END_2D 703 751 704 DO_2D( 0, 0, 0, 0 ) 752 DO_2D( 0, 0, 0, 0 ) ! no vector loop 705 753 706 754 ! tension**2 at T points … … 724 772 725 773 ! delta at T points 726 zdelta 727 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta) ) ! 0 if delta=0728 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch774 zdelta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 775 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0 776 pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch 729 777 730 778 END_2D … … 749 797 ! 5) diagnostics 750 798 !------------------------------------------------------------------------------! 751 DO_2D( 1, 1, 1, 1 )752 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice753 END_2D754 755 799 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 756 800 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 816 860 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 817 861 ENDIF 818 862 819 863 ! --- SIMIP --- ! 820 864 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 881 925 ENDIF 882 926 ! 927 ! --- convergence tests --- ! 928 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 929 IF( iom_use('uice_cvg') ) THEN 930 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 931 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 932 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 933 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 934 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 935 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 936 ENDIF 937 ENDIF 938 ENDIF 939 ! 940 DEALLOCATE( zmsk00, zmsk15 ) 941 ! 883 942 END SUBROUTINE ice_dyn_rhg_evp 943 944 945 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 946 !!---------------------------------------------------------------------- 947 !! *** ROUTINE rhg_cvg *** 948 !! 949 !! ** Purpose : check convergence of oce rheology 950 !! 951 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 952 !! during the sub timestepping of rheology so as: 953 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 954 !! This routine is called every sub-iteration, so it is cpu expensive 955 !! 956 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 957 !!---------------------------------------------------------------------- 958 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 959 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 960 !! 961 INTEGER :: it, idtime, istatus 962 INTEGER :: ji, jj ! dummy loop indices 963 REAL(wp) :: zresm ! local real 964 CHARACTER(len=20) :: clname 965 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 966 !!---------------------------------------------------------------------- 967 968 ! create file 969 IF( kt == nit000 .AND. kiter == 1 ) THEN 970 ! 971 IF( lwp ) THEN 972 WRITE(numout,*) 973 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 974 WRITE(numout,*) '~~~~~~~' 975 ENDIF 976 ! 977 IF( lwm ) THEN 978 clname = 'ice_cvg.nc' 979 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 980 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 981 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 982 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 983 istatus = NF90_ENDDEF(ncvgid) 984 ENDIF 985 ! 986 ENDIF 987 988 ! time 989 it = ( kt - 1 ) * kitermax + kiter 990 991 ! convergence 992 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 993 zresm = 0._wp 994 ELSE 995 DO_2D( 1, 1, 1, 1 ) 996 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 997 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 998 END_2D 999 zresm = MAXVAL( zres ) 1000 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1001 ENDIF 1002 1003 IF( lwm ) THEN 1004 ! write variables 1005 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1006 ! close file 1007 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 1008 ENDIF 1009 1010 END SUBROUTINE rhg_cvg 884 1011 885 1012 … … 939 1066 END SUBROUTINE rhg_evp_rst 940 1067 1068 941 1069 #else 942 1070 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
