- Timestamp:
- 2020-09-29T12:41:06+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13507 sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/ICE/icedyn_rhg_evp.F90
r12511 r13540 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 "do_loop_substitute.h90" 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 51 58 !!---------------------------------------------------------------------- 52 59 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 120 127 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 121 128 REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 129 REAl(wp) :: zbetau, zbetav 122 130 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 123 131 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars … … 126 134 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 127 135 ! 128 REAL(wp) :: zresm ! Maximal error on ice velocity129 136 REAL(wp) :: zintb, zintn ! dummy argument 130 137 REAL(wp) :: zfac_x, zfac_y … … 142 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 143 150 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 144 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence145 151 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 146 152 ! ! ocean surface (ssh_m) if ice is not embedded … … 156 162 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 157 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 158 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 159 165 160 166 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 161 167 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 162 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 163 171 !! --- diags 164 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00165 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 166 173 !! --- SIMIP diags … … 175 182 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 176 183 ! 177 !!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.... 178 192 !------------------------------------------------------------------------------! 179 193 ! 0) mask at F points for the ice 180 194 !------------------------------------------------------------------------------! 181 195 ! ocean/land mask 182 DO_2D _10_10196 DO_2D( 1, 0, 1, 0 ) 183 197 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 184 198 END_2D … … 186 200 187 201 ! Lateral boundary conditions on velocity (modify zfmask) 188 zwf(:,:) = zfmask(:,:) 189 DO_2D_00_00 202 DO_2D( 0, 0, 0, 0 ) 190 203 IF( zfmask(ji,jj) == 0._wp ) THEN 191 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) ) ) 204 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 205 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 192 206 ENDIF 193 207 END_2D 194 208 DO jj = 2, jpjm1 195 209 IF( zfmask(1,jj) == 0._wp ) THEN 196 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )210 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 197 211 ENDIF 198 212 IF( zfmask(jpi,jj) == 0._wp ) THEN 199 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )200 213 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 214 ENDIF 201 215 END DO 202 216 DO ji = 2, jpim1 203 217 IF( zfmask(ji,1) == 0._wp ) THEN 204 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )218 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 205 219 ENDIF 206 220 IF( zfmask(ji,jpj) == 0._wp ) THEN 207 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )221 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 208 222 ENDIF 209 223 END DO … … 219 233 z1_ecc2 = 1._wp / ecc2 220 234 221 ! Time step for subcycling222 zdtevp = rDt_ice / REAL( nn_nevp )223 z1_dtevp = 1._wp / zdtevp224 225 235 ! alpha parameters (Bouillon 2009) 226 236 IF( .NOT. ln_aEVP ) THEN 227 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 ) 228 239 zalph2 = zalph1 * z1_ecc2 229 240 230 241 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 231 242 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 243 ELSE 244 zdtevp = rdt_ice 245 ! zalpha parameters set later on adaptatively 232 246 ENDIF 247 z1_dtevp = 1._wp / zdtevp 233 248 234 249 ! Initialise stress tensor … … 241 256 242 257 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 243 IF( ln_landfast_L16 ) THEN ; zkt = rn_ tensile258 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 244 259 ELSE ; zkt = 0._wp 245 260 ENDIF … … 253 268 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 254 269 255 DO_2D _00_00270 DO_2D( 0, 0, 0, 0 ) 256 271 257 272 ! ice fraction at U-V points … … 299 314 300 315 END_2D 301 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1. , zdt_m, 'T', 1.)316 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp ) 302 317 ! 303 318 ! !== Landfast ice parameterization ==! 304 319 ! 305 320 IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 306 DO_2D _00_00321 DO_2D( 0, 0, 0, 0 ) 307 322 ! ice thickness at U-V points 308 323 zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 309 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) 310 325 ! ice-bottom stress at U points 311 zvCr = zaU(ji,jj) * rn_ depfra * hu(ji,jj,Kmm)312 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) ) ) 313 328 ! ice-bottom stress at V points 314 zvCr = zaV(ji,jj) * rn_ depfra * hv(ji,jj,Kmm)315 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) ) ) 316 331 ! ice_bottom stress at T points 317 zvCr = at_i(ji,jj) * rn_ depfra * ht(ji,jj)318 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) ) ) 319 334 END_2D 320 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )335 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp ) 321 336 ! 322 337 ELSE !-- no landfast 323 DO_2D _00_00338 DO_2D( 0, 0, 0, 0 ) 324 339 ztaux_base(ji,jj) = 0._wp 325 340 ztauy_base(ji,jj) = 0._wp … … 336 351 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 337 352 ! 338 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 339 !!$ DO jj = 1, jpjm1 340 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 341 !!$ zv_ice(:,jj) = v_ice(:,jj) 342 !!$ END DO 343 !!$ 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 344 360 345 361 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 346 DO_2D _10_10362 DO_2D( 1, 0, 1, 0 ) 347 363 348 364 ! shear at F points … … 352 368 353 369 END_2D 354 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )355 356 DO_2D _01_01370 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 371 372 DO_2D( 0, 1, 0, 1 ) ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 ! no vector loop 357 373 358 374 ! shear**2 at T points (doc eq. A16) … … 379 395 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 380 396 381 ! alpha & betafor aEVP397 ! alpha for aEVP 382 398 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 383 399 ! alpha = beta = sqrt(4*gamma) … … 387 403 zalph2 = zalph1 388 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 389 410 ENDIF 390 411 … … 394 415 395 416 END_2D 396 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 397 398 DO_2D_10_10 399 400 ! alpha & beta for aEVP 417 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 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 426 DO_2D( 1, 0, 1, 0 ) 427 428 ! alpha for aEVP 401 429 IF( ln_aEVP ) THEN 402 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) ) 403 431 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 404 zbeta(ji,jj) = zalph2 432 ! explicit: 433 ! z1_alph2 = 1._wp / zalph2 434 ! zalph2 = zalph2 - 1._wp 405 435 ENDIF 406 436 … … 414 444 415 445 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 416 DO_2D _00_00446 DO_2D( 0, 0, 0, 0 ) 417 447 ! !--- U points 418 448 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & … … 442 472 IF( MOD(jter,2) == 0 ) THEN ! even iterations 443 473 ! 444 DO_2D _00_00474 DO_2D( 0, 0, 0, 0 ) 445 475 ! !--- tau_io/(v_oce - v_ice) 446 476 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & … … 468 498 ! 469 499 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 470 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 471 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 472 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 473 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 474 & ) * 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 475 508 & ) * zmsk00y(ji,jj) 476 509 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 477 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity478 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)479 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast480 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0481 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin482 & ) 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) 483 516 ENDIF 484 517 END_2D 485 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )518 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 486 519 ! 487 520 #if defined key_agrif … … 491 524 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 492 525 ! 493 DO_2D _00_00526 DO_2D( 0, 0, 0, 0 ) 494 527 ! !--- tau_io/(u_oce - u_ice) 495 528 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & … … 517 550 ! 518 551 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 519 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 520 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 521 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 522 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 523 & ) * 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 524 560 & ) * zmsk00x(ji,jj) 525 561 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 526 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity527 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)528 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast529 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0530 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin531 & 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) 532 568 ENDIF 533 569 END_2D 534 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )570 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 535 571 ! 536 572 #if defined key_agrif … … 542 578 ELSE ! odd iterations 543 579 ! 544 DO_2D _00_00580 DO_2D( 0, 0, 0, 0 ) 545 581 ! !--- tau_io/(u_oce - u_ice) 546 582 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & … … 568 604 ! 569 605 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 570 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 571 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 572 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 573 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 574 & ) * 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 575 614 & ) * zmsk00x(ji,jj) 576 615 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 577 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj)* u_ice(ji,jj) & ! previous velocity578 & + zRHS + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)579 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast580 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0581 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin582 & 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) 583 622 ENDIF 584 623 END_2D 585 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )624 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 586 625 ! 587 626 #if defined key_agrif … … 591 630 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 592 631 ! 593 DO_2D _00_00632 DO_2D( 0, 0, 0, 0 ) 594 633 ! !--- tau_io/(v_oce - v_ice) 595 634 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & … … 617 656 ! 618 657 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 619 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 620 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 621 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 622 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 623 & ) * 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 624 666 & ) * zmsk00y(ji,jj) 625 667 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 626 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj)* v_ice(ji,jj) & ! previous velocity627 & + zRHS + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)628 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast629 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0630 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )& ! v_ice = v_oce/100 if mass < zmmin & conc < zamin631 & 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) 632 674 ENDIF 633 675 END_2D 634 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )676 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 635 677 ! 636 678 #if defined key_agrif … … 642 684 ENDIF 643 685 644 !!$ IF(sn_cfctl%l_prtctl) THEN ! Convergence test 645 !!$ DO jj = 2 , jpjm1 646 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 647 !!$ END DO 648 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 649 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 650 !!$ 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 ) 651 688 ! 652 689 ! ! ==================== ! 653 690 END DO ! end loop over jter ! 654 691 ! ! ==================== ! 692 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 655 693 ! 656 694 !------------------------------------------------------------------------------! 657 695 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 658 696 !------------------------------------------------------------------------------! 659 DO_2D _10_10697 DO_2D( 1, 0, 1, 0 ) 660 698 661 699 ! shear at F points … … 666 704 END_2D 667 705 668 DO_2D _00_00706 DO_2D( 0, 0, 0, 0 ) ! no vector loop 669 707 670 708 ! tension**2 at T points … … 693 731 694 732 END_2D 695 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1. , pdivu_i, 'T', 1., pdelta_i, 'T', 1.)733 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 696 734 697 735 ! --- Store the stress tensor for the next time step --- ! 698 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1. , zs2, 'T', 1., zs12, 'F', 1.)736 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp ) 699 737 pstress1_i (:,:) = zs1 (:,:) 700 738 pstress2_i (:,:) = zs2 (:,:) … … 705 743 ! 5) diagnostics 706 744 !------------------------------------------------------------------------------! 707 DO_2D_11_11708 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice709 END_2D710 711 745 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 712 746 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 713 747 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 714 748 ! 715 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1. , ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &716 & ztaux_bi, 'U', -1. , ztauy_bi, 'V', -1.)749 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, & 750 & ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 717 751 ! 718 752 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 734 768 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 735 769 ! 736 DO_2D _00_00770 DO_2D( 0, 0, 0, 0 ) 737 771 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 738 772 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & … … 751 785 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 752 786 END_2D 753 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1. , zsig2, 'T', 1., zsig3, 'T', 1.)787 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 754 788 ! 755 789 CALL iom_put( 'isig1' , zsig1 ) … … 763 797 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 764 798 ENDIF 765 799 766 800 ! --- SIMIP --- ! 767 801 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 768 802 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 769 803 ! 770 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1. , zspgV, 'V', -1., &771 & zCorU, 'U', -1. , zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1.)804 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, & 805 & zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 772 806 773 807 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) … … 785 819 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 786 820 ! 787 DO_2D _00_00821 DO_2D( 0, 0, 0, 0 ) 788 822 ! 2D ice mass, snow mass, area transport arrays (X, Y) 789 823 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) … … 801 835 END_2D 802 836 803 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1. , zdiag_ymtrp_ice, 'V', -1., &804 & zdiag_xmtrp_snw, 'U', -1. , zdiag_ymtrp_snw, 'V', -1., &805 & zdiag_xatrp , 'U', -1. , zdiag_yatrp , 'V', -1.)837 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, & 838 & zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, & 839 & zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp ) 806 840 807 841 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) … … 817 851 ENDIF 818 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 ! 819 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 820 937 821 938 … … 844 961 ! 845 962 IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist 846 CALL iom_get( numrir, jpdom_auto glo, 'stress1_i' , stress1_i)847 CALL iom_get( numrir, jpdom_auto glo, 'stress2_i' , stress2_i)848 CALL iom_get( numrir, jpdom_auto glo, 'stress12_i', stress12_i)963 CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' ) 964 CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' ) 965 CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' ) 849 966 ELSE ! start rheology from rest 850 967 IF(lwp) WRITE(numout,*) … … 875 992 END SUBROUTINE rhg_evp_rst 876 993 994 877 995 #else 878 996 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
