Changeset 12065 for NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2019-12-05T12:06:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icedyn_rhg_evp.F90
r10555 r12065 112 112 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 113 113 !! 114 LOGICAL, PARAMETER :: ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T)115 ! or only at the main time step (F)116 114 INTEGER :: ji, jj ! dummy loop indices 117 115 INTEGER :: jter ! local integers … … 123 121 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 124 122 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 125 REAL(wp) :: zTauO, zTauB, z TauE, zvel! temporary scalars123 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 126 124 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 127 125 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast … … 132 130 REAL(wp) :: zshear, zdum1, zdum2 133 131 ! 134 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors135 132 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 136 133 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 137 134 ! 138 135 REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points 139 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV! ice fraction on U/V points136 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 140 137 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 141 138 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 142 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points143 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ib , ztauV_ib ! ice-bottom stress at U-V points (landfast param)144 REAL(wp), DIMENSION(jpi,jpj) :: zspgU , zspgV ! surface pressure gradient at U/V points145 139 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 146 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses147 140 ! 148 141 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear … … 152 145 ! ! ocean surface (ssh_m) if ice is not embedded 153 146 ! ! ice bottom surface if ice is embedded 154 REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array 155 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array 156 ! 157 REAL(wp), DIMENSION(jpi,jpj) :: zswitchU, zswitchV ! dummy arrays 158 REAL(wp), DIMENSION(jpi,jpj) :: zmaskU, zmaskV ! mask for ice presence 147 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 148 REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points 149 REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array 150 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points 151 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points 152 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 153 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 154 ! 155 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 156 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 159 157 REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice 160 158 161 159 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 162 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 160 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 161 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 163 162 !! --- diags 164 REAL(wp), DIMENSION(jpi,jpj) :: z swi163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 165 164 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 166 165 !! --- SIMIP diags 167 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dx ! X-direction sea-surface tilt term (N/m2)170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dy ! X-direction sea-surface tilt term (N/m2)171 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstrx ! X-direction coriolis stress (N/m2)172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstry ! Y-direction coriolis stress (N/m2)173 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstrx ! X-direction internal stress (N/m2)174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstry ! Y-direction internal stress (N/m2)175 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_utau_oi ! X-direction ocean-ice stress176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress177 166 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 178 167 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) … … 254 243 CALL ice_strength 255 244 256 ! scale factors257 DO jj = 2, jpjm1258 DO ji = fs_2, fs_jpim1259 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )260 z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )261 END DO262 END DO263 264 245 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 265 IF( ln_landfast_L16 .OR. ln_landfast_home) THEN ; zkt = rn_tensile266 ELSE 246 IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile 247 ELSE ; zkt = 0._wp 267 248 ENDIF 268 249 ! … … 290 271 291 272 ! Ocean currents at U-V points 292 v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & 293 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 294 295 u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & 296 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 273 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) 274 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) 297 275 298 276 ! Coriolis at T points (m*f) … … 307 285 308 286 ! Drag ice-atm. 309 z TauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)310 z TauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)287 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 288 ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 311 289 312 290 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points … … 315 293 316 294 ! masks 317 zm askU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice318 zm askV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice295 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 296 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 319 297 320 298 ! switches 321 zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 322 zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 299 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 300 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 301 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 302 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 323 303 324 304 END DO … … 336 316 ! ice-bottom stress at U points 337 317 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 338 z TauU_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )318 ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 339 319 ! ice-bottom stress at V points 340 320 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 341 z TauV_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )321 ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 342 322 ! ice_bottom stress at T points 343 323 zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 344 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(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) ) ) 345 325 END DO 346 326 END DO 347 327 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 348 328 ! 349 ELSE IF( ln_landfast_home ) THEN !-- Home made329 ELSE !-- no landfast 350 330 DO jj = 2, jpjm1 351 331 DO ji = fs_2, fs_jpim1 352 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 353 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 354 END DO 355 END DO 356 ! 357 ELSE !-- no landfast 358 DO jj = 2, jpjm1 359 DO ji = fs_2, fs_jpim1 360 zTauU_ib(ji,jj) = 0._wp 361 zTauV_ib(ji,jj) = 0._wp 332 ztaux_base(ji,jj) = 0._wp 333 ztauy_base(ji,jj) = 0._wp 362 334 END DO 363 335 END DO 364 336 ENDIF 365 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )366 337 367 338 !------------------------------------------------------------------------------! … … 369 340 !------------------------------------------------------------------------------! 370 341 ! 371 ! ! ----------------------!342 ! ! ==================== ! 372 343 DO jter = 1 , nn_nevp ! loop over jter ! 373 ! ! ----------------------!344 ! ! ==================== ! 374 345 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 375 346 ! … … 476 447 & ) * r1_e1e2v(ji,jj) 477 448 ! 478 ! !--- u_ice atV point479 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &480 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)449 ! !--- ice currents at U-V point 450 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) 451 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) 481 452 ! 482 ! !--- v_ice at U point483 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &484 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)485 453 END DO 486 454 END DO … … 501 469 ! !--- tau_bottom/v_ice 502 470 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 503 zTauB = - zTauV_ib(ji,jj) / zvel 471 zTauB = ztauy_base(ji,jj) / zvel 472 ! !--- OceanBottom-to-Ice stress 473 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 504 474 ! 505 475 ! !--- Coriolis at V-points (energy conserving formulation) 506 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &476 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 507 477 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 508 478 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 509 479 ! 510 480 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 511 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 512 ! 513 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 514 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 481 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 482 ! 483 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 484 ! 1 = sliding friction : TauB < RHS 485 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 515 486 ! 516 487 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 517 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity518 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)519 & )/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast520 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0521 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin522 & ) * zmaskV(ji,jj)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 493 & ) * zmsk00y(ji,jj) 523 494 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 524 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity525 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)526 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast527 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0528 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin529 & ) * zmaskV(ji,jj)495 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 496 & + 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) + landfast 498 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 499 & ) * 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 & ) * zmsk00y(ji,jj) 530 501 ENDIF 531 502 END DO … … 537 508 CALL agrif_interp_ice( 'V' ) 538 509 #endif 539 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )510 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 540 511 ! 541 512 DO jj = 2, jpjm1 … … 549 520 ! !--- tau_bottom/u_ice 550 521 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 551 zTauB = - zTauU_ib(ji,jj) / zvel 522 zTauB = ztaux_base(ji,jj) / zvel 523 ! !--- OceanBottom-to-Ice stress 524 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 552 525 ! 553 526 ! !--- Coriolis at U-points (energy conserving formulation) 554 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &527 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 555 528 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 556 529 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 557 530 ! 558 531 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 559 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 560 ! 561 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 562 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 532 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 533 ! 534 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 535 ! 1 = sliding friction : TauB < RHS 536 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 563 537 ! 564 538 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 565 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity566 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)567 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast568 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0569 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin570 & ) * zmaskU(ji,jj)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 544 & ) * zmsk00x(ji,jj) 571 545 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 572 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity573 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)574 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast575 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0576 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin577 & ) * zmaskU(ji,jj)546 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 547 & + 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) + landfast 549 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 550 & ) * 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 551 & ) * zmsk00x(ji,jj) 578 552 ENDIF 579 553 END DO … … 585 559 CALL agrif_interp_ice( 'U' ) 586 560 #endif 587 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )561 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 588 562 ! 589 563 ELSE ! odd iterations … … 599 573 ! !--- tau_bottom/u_ice 600 574 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 601 zTauB = - zTauU_ib(ji,jj) / zvel 575 zTauB = ztaux_base(ji,jj) / zvel 576 ! !--- OceanBottom-to-Ice stress 577 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 602 578 ! 603 579 ! !--- Coriolis at U-points (energy conserving formulation) 604 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &580 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 605 581 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 606 582 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 607 583 ! 608 584 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 609 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 610 ! 611 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 612 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 585 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 586 ! 587 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 588 ! 1 = sliding friction : TauB < RHS 589 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 613 590 ! 614 591 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 615 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity616 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)617 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast618 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0619 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin620 & ) * zmaskU(ji,jj)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 597 & ) * zmsk00x(ji,jj) 621 598 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 622 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity623 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)624 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast625 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0626 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin627 & ) * zmaskU(ji,jj)599 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 600 & + 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) + landfast 602 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 603 & ) * 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 604 & ) * zmsk00x(ji,jj) 628 605 ENDIF 629 606 END DO … … 635 612 CALL agrif_interp_ice( 'U' ) 636 613 #endif 637 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )614 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 638 615 ! 639 616 DO jj = 2, jpjm1 … … 647 624 ! !--- tau_bottom/v_ice 648 625 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 649 zTauB = - zTauV_ib(ji,jj) / zvel 626 zTauB = ztauy_base(ji,jj) / zvel 627 ! !--- OceanBottom-to-Ice stress 628 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 650 629 ! 651 630 ! !--- Coriolis at v-points (energy conserving formulation) 652 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &631 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 653 632 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 654 633 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 655 634 ! 656 635 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 657 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 658 ! 659 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 660 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 636 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 637 ! 638 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 639 ! 1 = sliding friction : TauB < RHS 640 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 661 641 ! 662 642 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 663 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity664 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)665 & )/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast666 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0667 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin668 & ) * zmaskV(ji,jj)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 648 & ) * zmsk00y(ji,jj) 669 649 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 670 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity671 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)672 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast673 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0674 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin675 & ) * zmaskV(ji,jj)650 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 651 & + 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) + landfast 653 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 654 & ) * 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 655 & ) * zmsk00y(ji,jj) 676 656 ENDIF 677 657 END DO … … 683 663 CALL agrif_interp_ice( 'V' ) 684 664 #endif 685 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )665 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 686 666 ! 687 667 ENDIF … … 698 678 END DO ! end loop over jter ! 699 679 ! ! ==================== ! 700 !701 IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN702 CALL bdy_ice_dyn( 'U' )703 CALL bdy_ice_dyn( 'V' )704 ENDIF705 680 ! 706 681 !------------------------------------------------------------------------------! … … 761 736 DO jj = 1, jpj 762 737 DO ji = 1, jpi 763 z swi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice738 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 764 739 END DO 765 740 END DO 766 741 742 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 743 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 744 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 745 ! 746 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 747 & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 748 ! 749 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 750 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 751 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 752 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 753 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 754 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 755 ENDIF 756 767 757 ! --- divergence, shear and strength --- ! 768 IF( iom_use('icediv') ) CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:)) ! divergence769 IF( iom_use('iceshe') ) CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:)) ! shear770 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Icestrength771 772 ! --- charge ellipse--- !773 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN758 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 759 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 760 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 761 762 ! --- stress tensor --- ! 763 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 774 764 ! 775 765 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) … … 777 767 DO jj = 2, jpjm1 778 768 DO ji = 2, jpim1 779 zdum1 = ( z swi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point780 & z swi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) &781 & / MAX( 1._wp, z swi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )769 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 770 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 771 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 782 772 783 773 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 784 774 785 zdum2 = z swi(ji,jj) / MAX( 1._wp, strength(ji,jj) )775 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 786 776 787 777 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) … … 796 786 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 797 787 ! 798 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) 799 IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) 800 IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) 801 ! 788 CALL iom_put( 'isig1' , zsig1 ) 789 CALL iom_put( 'isig2' , zsig2 ) 790 CALL iom_put( 'isig3' , zsig3 ) 791 ! 792 ! Stress tensor invariants (normal and shear stress N/m) 793 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 794 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 795 802 796 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 803 797 ENDIF 804 798 805 799 ! --- SIMIP --- ! 806 IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & 807 & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & 808 & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 809 & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN 810 811 ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & 812 & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & 813 & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 814 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) 815 800 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 801 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 802 ! 803 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 804 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 805 806 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 807 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 808 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 809 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 810 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 811 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 812 ENDIF 813 814 IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 815 & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 816 ! 817 ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 818 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 819 ! 816 820 DO jj = 2, jpjm1 817 821 DO ji = 2, jpim1 818 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice819 820 ! Stress tensor invariants (normal and shear stress N/m)821 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress822 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress823 824 ! Stress terms of the momentum equation (N/m2)825 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term826 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch827 828 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term829 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch830 831 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term832 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch833 834 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress835 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch836 837 822 ! 2D ice mass, snow mass, area transport arrays (X, Y) 838 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch839 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch840 823 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 824 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 825 841 826 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 842 827 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 843 828 844 829 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 845 830 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 846 831 847 832 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 848 833 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 849 850 END DO 851 END DO 852 853 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 854 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 855 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 856 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 857 858 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi , 'U', -1., zdiag_vtau_oi , 'V', -1., & 859 & zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 860 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 861 & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) 862 863 IF( iom_use('normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress 864 IF( iom_use('sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress 865 IF( iom_use('dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) 866 IF( iom_use('dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) 867 IF( iom_use('corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) 868 IF( iom_use('corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) 869 IF( iom_use('intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) 870 IF( iom_use('intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) 871 IF( iom_use('utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) 872 IF( iom_use('vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) 873 IF( iom_use('xmtrpice') ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) 874 IF( iom_use('ymtrpice') ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport 875 IF( iom_use('xmtrpsnw') ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) 876 IF( iom_use('ymtrpsnw') ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport 877 IF( iom_use('xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport 878 IF( iom_use('yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport 879 880 DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & 881 & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & 882 & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 883 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 834 835 END DO 836 END DO 837 838 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 839 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 840 & zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. ) 841 842 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) 843 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 844 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 845 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport 846 CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport 847 CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport 848 849 DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 850 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 884 851 885 852 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
