- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ICE/icedyn_rhg_evp.F90
r14433 r15548 48 48 PUBLIC rhg_evp_rst ! called by icedyn_rhg.F90 49 49 50 !! for convergence tests 51 INTEGER :: ncvgid ! netcdf file id 52 INTEGER :: nvarid ! netcdf variable id 53 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice 54 50 55 !! * Substitutions 51 56 # include "do_loop_substitute.h90" 52 57 # include "domzgr_substitute.h90" 53 54 !! for convergence tests55 INTEGER :: ncvgid ! netcdf file id56 INTEGER :: nvarid ! netcdf variable id57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk1558 58 !!---------------------------------------------------------------------- 59 59 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 134 134 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 135 135 ! 136 REAL(wp) :: zintb, zintn ! dummy argument137 136 REAL(wp) :: zfac_x, zfac_y 138 REAL(wp) :: zshear, zdum1, zdum2139 137 ! 140 138 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points … … 161 159 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 162 160 ! 161 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15 163 162 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 164 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 165 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice166 164 167 165 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter … … 180 178 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 181 179 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 180 !! -- advect fields at the rheology time step for the calculation of strength 181 !! it seems that convergence is worse when ll_advups=true. So it not really a good idea 182 LOGICAL :: ll_advups = .FALSE. 183 REAL(wp) :: zdt_ups 184 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmp 185 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: za_i_ups, zv_i_ups ! tracers advected upstream 182 186 !!------------------------------------------------------------------- 183 187 … … 185 189 ! 186 190 ! for diagnostics and convergence tests 187 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 188 DO_2D( 1, 1, 1, 1 ) 191 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 189 192 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 190 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less191 193 END_2D 192 ! 193 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 194 IF( nn_rhg_chkcvg > 0 ) THEN 195 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 196 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 197 END_2D 198 ENDIF 199 ! 194 200 !------------------------------------------------------------------------------! 195 201 ! 0) mask at F points for the ice 196 202 !------------------------------------------------------------------------------! 197 ! ocean/land mask 198 DO_2D( 1, 0, 1, 0 ) 199 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 200 END_2D 201 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp) 202 203 ! Lateral boundary conditions on velocity (modify zfmask) 204 DO_2D( 0, 0, 0, 0 ) 205 IF( zfmask(ji,jj) == 0._wp ) THEN 206 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 207 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 203 IF( kt == nit000 ) THEN 204 ! ocean/land mask 205 ALLOCATE( fimask(jpi,jpj) ) 206 IF( rn_ishlat == 0._wp ) THEN 207 DO_2D( 0, 0, 0, 0 ) 208 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 209 END_2D 210 ELSE 211 DO_2D( 0, 0, 0, 0 ) 212 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 213 ! Lateral boundary conditions on velocity (modify fimask) 214 IF( fimask(ji,jj) == 0._wp ) THEN 215 fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 216 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 217 ENDIF 218 END_2D 208 219 ENDIF 209 END_2D 210 DO jj = 2, jpjm1 211 IF( zfmask(1,jj) == 0._wp ) THEN 212 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 213 ENDIF 214 IF( zfmask(jpi,jj) == 0._wp ) THEN 215 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 216 ENDIF 217 END DO 218 DO ji = 2, jpim1 219 IF( zfmask(ji,1) == 0._wp ) THEN 220 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 221 ENDIF 222 IF( zfmask(ji,jpj) == 0._wp ) THEN 223 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 224 ENDIF 225 END DO 226 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) 227 220 CALL lbc_lnk( 'icedyn_rhg_evp', fimask, 'F', 1._wp ) 221 ENDIF 228 222 !------------------------------------------------------------------------------! 229 223 ! 1) define some variables and initialize arrays … … 244 238 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 245 239 ELSE 246 zdtevp = r dt_ice240 zdtevp = rDt_ice 247 241 ! zalpha parameters set later on adaptatively 248 242 ENDIF … … 270 264 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 271 265 272 DO_2D( 0, 0, 0, 0 ) 266 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 267 zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ) ! Ice/snow mass at U-V points 268 zmf (ji,jj) = zm1 * ff_t(ji,jj) ! Coriolis at T points (m*f) 269 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) ! dt/m at T points (for alpha and beta coefficients) 270 END_2D 271 272 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 273 273 274 274 ! ice fraction at U-V points … … 284 284 285 285 ! Ocean currents at U-V points 286 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 287 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 288 289 ! Coriolis at T points (m*f) 290 zmf(ji,jj) = zm1 * ff_t(ji,jj) 291 292 ! dt/m at T points (for alpha and beta coefficients) 293 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 286 ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) 287 v_oceU(ji,jj) = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1) 288 u_oceV(ji,jj) = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1) 294 289 295 290 ! m/dt … … 316 311 317 312 END_2D 318 CALL lbc_lnk( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )319 313 ! 320 314 ! !== Landfast ice parameterization ==! 321 315 ! 322 316 IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 323 DO_2D( 0, 0, 0, 0)317 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 324 318 ! ice thickness at U-V points 325 319 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) … … 338 332 ! 339 333 ELSE !-- no landfast 340 DO_2D( 0, 0, 0, 0)334 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 341 335 ztaux_base(ji,jj) = 0._wp 342 336 ztauy_base(ji,jj) = 0._wp … … 362 356 363 357 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 364 DO_2D( 1, 0, 1, 0)358 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 365 359 366 360 ! shear at F points 367 361 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 368 362 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 369 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)363 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 370 364 371 365 END_2D … … 393 387 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 394 388 389 ! P/delta at T points 390 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 391 395 392 END_2D 396 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 397 398 ! P/delta at T points 399 DO_2D( 1, 1, 1, 1 ) 400 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 401 END_2D 402 403 DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 393 CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp ) 394 395 ! 396 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 404 397 405 398 ! divergence at T points (duplication to avoid communications) 406 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 407 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 399 ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) 400 zdiv = ( (e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)) & 401 & + (e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)) & 408 402 & ) * r1_e1e2t(ji,jj) 409 403 … … 436 430 ! Save beta at T-points for further computations 437 431 IF( ln_aEVP ) THEN 438 DO_2D( 1, 1, 1, 1)432 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 439 433 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 440 434 END_2D 441 435 ENDIF 442 436 443 DO_2D( 1, 0, 1, 0)437 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 444 438 445 439 ! alpha for aEVP … … 453 447 454 448 ! P/delta at F points 455 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 449 ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) 450 zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) ) 456 451 457 452 ! stress at F points (zkt/=0 if landfast) … … 461 456 462 457 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 463 DO_2D( 0, 0, 0, 0 ) 458 ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) 459 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 464 460 ! !--- U points 465 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &461 zfU(ji,jj) = 0.5_wp * ( (( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 466 462 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 467 & ) * r1_e2u(ji,jj) &463 & ) * r1_e2u(ji,jj)) & 468 464 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 469 465 & ) * 2._wp * r1_e1u(ji,jj) & … … 471 467 ! 472 468 ! !--- V points 473 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &469 zfV(ji,jj) = 0.5_wp * ( (( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 474 470 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 475 & ) * r1_e1v(ji,jj) &471 & ) * r1_e1v(ji,jj)) & 476 472 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 477 473 & ) * 2._wp * r1_e2v(ji,jj) & … … 479 475 ! 480 476 ! !--- ice currents at U-V point 481 v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)482 u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)477 v_iceU(ji,jj) = 0.25_wp * ( (v_ice(ji,jj) + v_ice(ji,jj-1)) + (v_ice(ji+1,jj) + v_ice(ji+1,jj-1)) ) * umask(ji,jj,1) 478 u_iceV(ji,jj) = 0.25_wp * ( (u_ice(ji,jj) + u_ice(ji-1,jj)) + (u_ice(ji,jj+1) + u_ice(ji-1,jj+1)) ) * vmask(ji,jj,1) 483 479 ! 484 480 END_2D … … 489 485 IF( MOD(jter,2) == 0 ) THEN ! even iterations 490 486 ! 491 DO_2D( 0, 0, 0, 0)487 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 492 488 ! !--- tau_io/(v_oce - v_ice) 493 489 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & … … 533 529 ENDIF 534 530 END_2D 535 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 536 ! 537 #if defined key_agrif 538 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 539 CALL agrif_interp_ice( 'V' ) 540 #endif 541 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 531 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 542 532 ! 543 533 DO_2D( 0, 0, 0, 0 ) … … 585 575 ENDIF 586 576 END_2D 587 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 588 ! 589 #if defined key_agrif 590 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 591 CALL agrif_interp_ice( 'U' ) 592 #endif 593 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 577 IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 578 ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 579 ENDIF 594 580 ! 595 581 ELSE ! odd iterations 596 582 ! 597 DO_2D( 0, 0, 0, 0)583 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 598 584 ! !--- tau_io/(u_oce - u_ice) 599 585 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & … … 639 625 ENDIF 640 626 END_2D 641 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 642 ! 643 #if defined key_agrif 644 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 645 CALL agrif_interp_ice( 'U' ) 646 #endif 647 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 627 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) 648 628 ! 649 629 DO_2D( 0, 0, 0, 0 ) … … 691 671 ENDIF 692 672 END_2D 693 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 694 ! 673 IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) 674 ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 675 ENDIF 676 ! 677 ENDIF 678 ! 695 679 #if defined key_agrif 696 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 697 CALL agrif_interp_ice( 'V' ) 680 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 681 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 682 CALL agrif_interp_ice( 'U' ) 683 CALL agrif_interp_ice( 'V' ) 698 684 #endif 699 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 700 ! 685 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 686 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 687 ! 688 ! convergence test 689 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) 690 ! 691 ! 692 ! --- change strength according to advected a_i and v_i (upstream for now) --- ! 693 IF( ll_advups .AND. ln_str_H79 ) THEN 694 ! 695 IF( jter == 1 ) THEN ! init 696 ALLOCATE( za_i_ups(jpi,jpj,jpl), zv_i_ups(jpi,jpj,jpl) ) 697 ALLOCATE( ztmp(jpi,jpj) ) 698 zdt_ups = rDt_ice / REAL( nn_nevp ) 699 za_i_ups(:,:,:) = a_i(:,:,:) 700 zv_i_ups(:,:,:) = v_i(:,:,:) 701 ELSE 702 CALL lbc_lnk( 'icedyn_rhg_evp', za_i_ups, 'T', 1.0_wp, zv_i_ups, 'T', 1.0_wp ) 703 ENDIF 704 ! 705 CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, za_i_ups ) ! upstream advection: a_i 706 CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, zv_i_ups ) ! upstream advection: v_i 707 ! 708 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! strength 709 strength(ji,jj) = rn_pstar * SUM( zv_i_ups(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - SUM( za_i_ups(ji,jj,:) ) ) ) 710 END_2D 711 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp ) 712 ! 713 DO_2D( 0, 0, 0, 0 ) ! strength smoothing 714 IF( SUM( za_i_ups(ji,jj,:) ) > 0._wp ) THEN 715 ztmp(ji,jj) = ( 4._wp * strength(ji,jj) + strength(ji-1,jj ) + strength(ji+1,jj ) & 716 & + strength(ji ,jj-1) + strength(ji ,jj+1) & 717 & ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 718 ELSE 719 ztmp(ji,jj) = 0._wp 720 ENDIF 721 END_2D 722 DO_2D( 0, 0, 0, 0 ) 723 strength(ji,jj) = ztmp(ji,jj) 724 END_2D 725 ! 726 IF( jter == nn_nevp ) THEN 727 DEALLOCATE( za_i_ups, zv_i_ups ) 728 DEALLOCATE( ztmp ) 729 ENDIF 701 730 ENDIF 702 703 ! convergence test704 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )705 !706 731 ! ! ==================== ! 707 732 END DO ! end loop over jter ! … … 709 734 IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) 710 735 ! 736 IF( ll_advups .AND. ln_str_H79 ) CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp ) 737 ! 711 738 !------------------------------------------------------------------------------! 712 739 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 713 740 !------------------------------------------------------------------------------! 714 DO_2D( 1, 0, 1, 0)741 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 715 742 716 743 ! shear at F points 717 744 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 718 745 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 719 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)746 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 720 747 721 748 END_2D … … 782 809 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 783 810 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 811 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 784 812 785 813 ! --- Stress tensor invariants (SIMIP diags) --- ! … … 788 816 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 789 817 ! 790 DO_2D( 1, 1, 1, 1)818 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 791 819 792 820 ! Ice stresses … … 800 828 801 829 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 802 zsig_I (ji,jj) = zsig1 * 0.5_wp 803 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress830 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 831 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 804 832 805 833 END_2D … … 821 849 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 822 850 ! 823 DO_2D( 1, 1, 1, 1)851 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 824 852 825 853 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates … … 832 860 833 861 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 834 zsig_I(ji,jj) = zsig1 * 0.5_wp 835 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress862 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 863 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 836 864 837 865 ! Normalized principal stresses (used to display the ellipse) … … 914 942 ENDIF 915 943 ! 916 DEALLOCATE( zmsk00, zmsk15 )917 !918 944 END SUBROUTINE ice_dyn_rhg_evp 919 945 920 946 921 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )947 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 ) 922 948 !!---------------------------------------------------------------------- 923 949 !! *** ROUTINE rhg_cvg *** … … 934 960 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 935 961 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 962 REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15 936 963 !! 937 964 INTEGER :: it, idtime, istatus … … 939 966 REAL(wp) :: zresm ! local real 940 967 CHARACTER(len=20) :: clname 941 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 968 LOGICAL :: ll_maxcvg 969 REAL(wp), DIMENSION(jpi,jpj,2) :: zres 970 REAL(wp), DIMENSION(2) :: ztmp 942 971 !!---------------------------------------------------------------------- 943 972 ll_maxcvg = .FALSE. 973 ! 944 974 ! create file 945 975 IF( kt == nit000 .AND. kiter == 1 ) THEN … … 956 986 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 957 987 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 958 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE 988 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 959 989 istatus = NF90_ENDDEF(ncvgid) 960 990 ENDIF … … 963 993 964 994 ! time 965 it = ( kt - 1) * kitermax + kiter995 it = ( kt - nit000 ) * kitermax + kiter 966 996 967 997 ! convergence … … 969 999 zresm = 0._wp 970 1000 ELSE 971 DO_2D( 1, 1, 1, 1 ) 972 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 973 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 974 END_2D 975 zresm = MAXVAL( zres ) 976 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1001 zresm = 0._wp 1002 IF( ll_maxcvg ) THEN ! error max over the domain 1003 DO_2D( 0, 0, 0, 0 ) 1004 zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1005 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 1006 END_2D 1007 CALL mpp_max( 'icedyn_rhg_evp', zresm ) 1008 ELSE ! error averaged over the domain 1009 DO_2D( 0, 0, 0, 0 ) 1010 zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1011 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) 1012 zres(ji,jj,2) = pmsk15(ji,jj) 1013 END_2D 1014 ztmp(:) = glob_sum_vec( 'icedyn_rhg_evp', zres ) 1015 IF( ztmp(2) /= 0._wp ) zresm = ztmp(1) / ztmp(2) 1016 ENDIF 977 1017 ENDIF 978 1018 … … 981 1021 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 982 1022 ! close file 983 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid)1023 IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) 984 1024 ENDIF 985 1025 … … 1042 1082 END SUBROUTINE rhg_evp_rst 1043 1083 1084 SUBROUTINE rhg_upstream( jter, pdt, pu, pv, pt ) 1085 !!--------------------------------------------------------------------- 1086 !! *** ROUTINE rhg_upstream *** 1087 !! 1088 !! ** Purpose : compute the upstream fluxes and upstream guess of tracer 1089 !!---------------------------------------------------------------------- 1090 INTEGER , INTENT(in ) :: jter 1091 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1092 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components 1093 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer fields 1094 ! 1095 INTEGER :: ji, jj, jl ! dummy loop indices 1096 REAL(wp) :: ztra ! local scalar 1097 LOGICAL :: ll_upsxy = .TRUE. 1098 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zpt ! upstream fluxes and tracer guess 1099 !!---------------------------------------------------------------------- 1100 DO jl = 1, jpl 1101 IF( .NOT. ll_upsxy ) THEN !** no alternate directions **! 1102 ! 1103 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 1104 zfu_ups(ji,jj) = MAX(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji+1,jj,jl) 1105 zfv_ups(ji,jj) = MAX(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj+1,jl) 1106 END_2D 1107 ! 1108 ELSE !** alternate directions **! 1109 ! 1110 IF( MOD(jter,2) == 1 ) THEN !== odd ice time step: adv_x then adv_y ==! 1111 ! 1112 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls ) !-- flux in x-direction 1113 zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji ,jj,jl) + & 1114 & MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 1115 END_2D 1116 ! 1117 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls ) !-- first guess of tracer from u-flux 1118 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) 1119 zpt(ji,jj) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1120 END_2D 1121 ! 1122 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) !-- flux in y-direction 1123 zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj ) + & 1124 & MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj+1) 1125 END_2D 1126 ! 1127 ELSE !== even ice time step: adv_y then adv_x ==! 1128 ! 1129 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 ) !-- flux in y-direction 1130 zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj ,jl) + & 1131 & MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 1132 END_2D 1133 ! 1134 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) !-- first guess of tracer from v-flux 1135 ztra = - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 1136 zpt(ji,jj) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1137 END_2D 1138 ! 1139 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) !-- flux in x-direction 1140 zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji ,jj) + & 1141 & MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji+1,jj) 1142 END_2D 1143 ! 1144 ENDIF 1145 ! 1146 ENDIF 1147 ! 1148 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 1149 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 1150 pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 1151 END_2D 1152 END DO 1153 ! 1154 END SUBROUTINE rhg_upstream 1044 1155 1045 1156 #else
Note: See TracChangeset
for help on using the changeset viewer.
