Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rhg_evp.F90
r11480 r11822 113 113 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 114 114 !! 115 LOGICAL, PARAMETER :: ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T)116 ! or only at the main time step (F)117 115 INTEGER :: ji, jj ! dummy loop indices 118 116 INTEGER :: jter ! local integers … … 124 122 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 125 123 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 126 REAL(wp) :: zTauO, zTauB, z TauE, zvel! temporary scalars124 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 127 125 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 128 126 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast … … 133 131 REAL(wp) :: zshear, zdum1, zdum2 134 132 ! 135 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors136 133 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 137 134 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 138 135 ! 139 136 REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points 140 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV! ice fraction on U/V points137 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 141 138 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 142 139 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 143 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points144 REAL(wp), DIMENSION(jpi,jpj) :: zTauU_ib , ztauV_ib ! ice-bottom stress at U-V points (landfast param)145 REAL(wp), DIMENSION(jpi,jpj) :: zspgU , zspgV ! surface pressure gradient at U/V points146 140 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 147 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses148 141 ! 149 142 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear … … 153 146 ! ! ocean surface (ssh_m) if ice is not embedded 154 147 ! ! ice bottom surface if ice is embedded 155 REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array 156 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array 157 ! 158 REAL(wp), DIMENSION(jpi,jpj) :: zswitchU, zswitchV ! dummy arrays 159 REAL(wp), DIMENSION(jpi,jpj) :: zmaskU, zmaskV ! mask for ice presence 148 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 149 REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points 150 REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array 151 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points 152 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points 153 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 154 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 155 ! 156 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 157 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 160 158 REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice 161 159 162 160 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 163 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 161 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 162 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 164 163 !! --- diags 165 REAL(wp), DIMENSION(jpi,jpj) :: z swi164 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 166 165 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 167 166 !! --- SIMIP diags 168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dx ! X-direction sea-surface tilt term (N/m2)171 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dy ! X-direction sea-surface tilt term (N/m2)172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstrx ! X-direction coriolis stress (N/m2)173 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstry ! Y-direction coriolis stress (N/m2)174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstrx ! X-direction internal stress (N/m2)175 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstry ! Y-direction internal stress (N/m2)176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_utau_oi ! X-direction ocean-ice stress177 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress178 167 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 179 168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) … … 255 244 CALL ice_strength 256 245 257 ! scale factors258 DO jj = 2, jpjm1259 DO ji = fs_2, fs_jpim1260 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )261 z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )262 END DO263 END DO264 265 246 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 266 IF( ln_landfast_L16 .OR. ln_landfast_home) THEN ; zkt = rn_tensile267 ELSE 247 IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile 248 ELSE ; zkt = 0._wp 268 249 ENDIF 269 250 ! … … 291 272 292 273 ! Ocean currents at U-V points 293 v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & 294 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 295 296 u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & 297 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 274 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) 275 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) 298 276 299 277 ! Coriolis at T points (m*f) … … 308 286 309 287 ! Drag ice-atm. 310 z TauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)311 z TauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)288 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 289 ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 312 290 313 291 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points … … 316 294 317 295 ! masks 318 zm askU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice319 zm askV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice296 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 297 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 320 298 321 299 ! switches 322 zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 323 zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 300 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 301 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 302 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 303 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 324 304 325 305 END DO … … 337 317 ! ice-bottom stress at U points 338 318 zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm) 339 z TauU_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )319 ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 340 320 ! ice-bottom stress at V points 341 321 zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm) 342 z TauV_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )322 ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 343 323 ! ice_bottom stress at T points 344 324 zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj) 345 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )325 tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 346 326 END DO 347 327 END DO 348 328 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 349 329 ! 350 ELSE IF( ln_landfast_home ) THEN !-- Home made330 ELSE !-- no landfast 351 331 DO jj = 2, jpjm1 352 332 DO ji = fs_2, fs_jpim1 353 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 354 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 355 END DO 356 END DO 357 ! 358 ELSE !-- no landfast 359 DO jj = 2, jpjm1 360 DO ji = fs_2, fs_jpim1 361 zTauU_ib(ji,jj) = 0._wp 362 zTauV_ib(ji,jj) = 0._wp 333 ztaux_base(ji,jj) = 0._wp 334 ztauy_base(ji,jj) = 0._wp 363 335 END DO 364 336 END DO 365 337 ENDIF 366 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )367 338 368 339 !------------------------------------------------------------------------------! … … 370 341 !------------------------------------------------------------------------------! 371 342 ! 372 ! ! ----------------------!343 ! ! ==================== ! 373 344 DO jter = 1 , nn_nevp ! loop over jter ! 374 ! ! ----------------------!345 ! ! ==================== ! 375 346 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 376 347 ! … … 477 448 & ) * r1_e1e2v(ji,jj) 478 449 ! 479 ! !--- u_ice atV point480 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &481 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)450 ! !--- ice currents at U-V point 451 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) 452 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) 482 453 ! 483 ! !--- v_ice at U point484 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &485 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)486 454 END DO 487 455 END DO … … 502 470 ! !--- tau_bottom/v_ice 503 471 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 504 zTauB = - zTauV_ib(ji,jj) / zvel 472 zTauB = ztauy_base(ji,jj) / zvel 473 ! !--- OceanBottom-to-Ice stress 474 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 505 475 ! 506 476 ! !--- Coriolis at V-points (energy conserving formulation) 507 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &477 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 508 478 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 509 479 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 510 480 ! 511 481 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 512 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 513 ! 514 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 515 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 482 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 483 ! 484 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 485 ! 1 = sliding friction : TauB < RHS 486 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 516 487 ! 517 488 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 518 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity519 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)520 & )/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast521 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0522 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin523 & ) * zmaskV(ji,jj)489 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 490 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 491 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 492 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 493 & ) * 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 494 & ) * zmsk00y(ji,jj) 524 495 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 525 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity526 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)527 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast528 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0529 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin530 & ) * zmaskV(ji,jj)496 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 497 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 498 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 499 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 500 & ) * 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 501 & ) * zmsk00y(ji,jj) 531 502 ENDIF 532 503 END DO … … 538 509 CALL agrif_interp_ice( 'V' ) 539 510 #endif 540 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )511 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 541 512 ! 542 513 DO jj = 2, jpjm1 … … 550 521 ! !--- tau_bottom/u_ice 551 522 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 552 zTauB = - zTauU_ib(ji,jj) / zvel 523 zTauB = ztaux_base(ji,jj) / zvel 524 ! !--- OceanBottom-to-Ice stress 525 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 553 526 ! 554 527 ! !--- Coriolis at U-points (energy conserving formulation) 555 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &528 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 556 529 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 557 530 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 558 531 ! 559 532 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 560 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 561 ! 562 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 563 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 533 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 534 ! 535 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 536 ! 1 = sliding friction : TauB < RHS 537 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 564 538 ! 565 539 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 566 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity567 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)568 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast569 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0570 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin571 & ) * zmaskU(ji,jj)540 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 541 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 542 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 543 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 544 & ) * 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 545 & ) * zmsk00x(ji,jj) 572 546 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 573 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity574 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)575 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast576 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0577 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin578 & ) * zmaskU(ji,jj)547 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 548 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 549 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 550 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 551 & ) * 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 & ) * zmsk00x(ji,jj) 579 553 ENDIF 580 554 END DO … … 586 560 CALL agrif_interp_ice( 'U' ) 587 561 #endif 588 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )562 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 589 563 ! 590 564 ELSE ! odd iterations … … 600 574 ! !--- tau_bottom/u_ice 601 575 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 602 zTauB = - zTauU_ib(ji,jj) / zvel 576 zTauB = ztaux_base(ji,jj) / zvel 577 ! !--- OceanBottom-to-Ice stress 578 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 603 579 ! 604 580 ! !--- Coriolis at U-points (energy conserving formulation) 605 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &581 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 606 582 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 607 583 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 608 584 ! 609 585 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 610 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 611 ! 612 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 613 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 586 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 587 ! 588 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 589 ! 1 = sliding friction : TauB < RHS 590 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 614 591 ! 615 592 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 616 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity617 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)618 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast619 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0620 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin621 & ) * zmaskU(ji,jj)593 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 594 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 595 & / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 596 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 597 & ) * 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 598 & ) * zmsk00x(ji,jj) 622 599 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 623 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity624 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)625 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast626 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0627 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin628 & ) * zmaskU(ji,jj)600 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 601 & + zRHS + zTauO * u_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 602 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 603 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 604 & ) * 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 605 & ) * zmsk00x(ji,jj) 629 606 ENDIF 630 607 END DO … … 636 613 CALL agrif_interp_ice( 'U' ) 637 614 #endif 638 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )615 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 639 616 ! 640 617 DO jj = 2, jpjm1 … … 648 625 ! !--- tau_bottom/v_ice 649 626 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 650 zTauB = - zTauV_ib(ji,jj) / zvel 627 zTauB = ztauy_base(ji,jj) / zvel 628 ! !--- OceanBottom-to-Ice stress 629 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 651 630 ! 652 631 ! !--- Coriolis at v-points (energy conserving formulation) 653 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &632 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 654 633 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 655 634 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 656 635 ! 657 636 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 658 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 659 ! 660 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 661 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 637 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 638 ! 639 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 640 ! 1 = sliding friction : TauB < RHS 641 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 662 642 ! 663 643 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 664 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity665 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)666 & )/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast667 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0668 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin669 & ) * zmaskV(ji,jj)644 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 645 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 646 & / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 647 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 648 & ) * 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 649 & ) * zmsk00y(ji,jj) 670 650 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 671 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity672 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)673 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast674 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0675 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin676 & ) * zmaskV(ji,jj)651 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 652 & + zRHS + zTauO * v_ice(ji,jj) ) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 653 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 654 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 655 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 656 & ) * zmsk00y(ji,jj) 677 657 ENDIF 678 658 END DO … … 684 664 CALL agrif_interp_ice( 'V' ) 685 665 #endif 686 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )666 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 687 667 ! 688 668 ENDIF … … 699 679 END DO ! end loop over jter ! 700 680 ! ! ==================== ! 701 !702 IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN703 CALL bdy_ice_dyn( 'U' )704 CALL bdy_ice_dyn( 'V' )705 ENDIF706 681 ! 707 682 !------------------------------------------------------------------------------! … … 762 737 DO jj = 1, jpj 763 738 DO ji = 1, jpi 764 z swi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice739 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 765 740 END DO 766 741 END DO 767 742 743 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 744 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 745 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 746 ! 747 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 748 & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 749 ! 750 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 751 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 752 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 753 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 754 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 755 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 756 ENDIF 757 768 758 ! --- divergence, shear and strength --- ! 769 IF( iom_use('icediv') ) CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:)) ! divergence770 IF( iom_use('iceshe') ) CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:)) ! shear771 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Icestrength772 773 ! --- charge ellipse--- !774 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN759 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 760 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 761 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 762 763 ! --- stress tensor --- ! 764 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 775 765 ! 776 766 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) … … 778 768 DO jj = 2, jpjm1 779 769 DO ji = 2, jpim1 780 zdum1 = ( z swi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point781 & z swi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) &782 & / MAX( 1._wp, z swi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )770 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 771 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 772 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 783 773 784 774 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 785 775 786 zdum2 = z swi(ji,jj) / MAX( 1._wp, strength(ji,jj) )776 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 787 777 788 778 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) … … 797 787 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 798 788 ! 799 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) 800 IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) 801 IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) 802 ! 789 CALL iom_put( 'isig1' , zsig1 ) 790 CALL iom_put( 'isig2' , zsig2 ) 791 CALL iom_put( 'isig3' , zsig3 ) 792 ! 793 ! Stress tensor invariants (normal and shear stress N/m) 794 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 795 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 796 803 797 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 804 798 ENDIF 805 799 806 800 ! --- SIMIP --- ! 807 IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & 808 & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & 809 & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 810 & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN 811 812 ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & 813 & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & 814 & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 815 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) 816 801 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 802 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 803 ! 804 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 805 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 806 807 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 808 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 809 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 810 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 811 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 812 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 813 ENDIF 814 815 IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 816 & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 817 ! 818 ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 819 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 820 ! 817 821 DO jj = 2, jpjm1 818 822 DO ji = 2, jpim1 819 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice820 821 ! Stress tensor invariants (normal and shear stress N/m)822 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress823 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress824 825 ! Stress terms of the momentum equation (N/m2)826 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term827 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch828 829 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term830 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch831 832 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term833 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch834 835 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress836 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch837 838 823 ! 2D ice mass, snow mass, area transport arrays (X, Y) 839 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch840 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch841 824 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 825 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 826 842 827 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 843 828 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 844 829 845 830 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 846 831 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 847 832 848 833 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 849 834 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 850 851 END DO 852 END DO 853 854 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 855 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 856 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 857 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 858 859 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi , 'U', -1., zdiag_vtau_oi , 'V', -1., & 860 & zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 861 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 862 & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) 863 864 IF( iom_use('normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress 865 IF( iom_use('sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress 866 IF( iom_use('dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) 867 IF( iom_use('dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) 868 IF( iom_use('corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) 869 IF( iom_use('corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) 870 IF( iom_use('intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) 871 IF( iom_use('intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) 872 IF( iom_use('utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) 873 IF( iom_use('vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) 874 IF( iom_use('xmtrpice') ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) 875 IF( iom_use('ymtrpice') ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport 876 IF( iom_use('xmtrpsnw') ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) 877 IF( iom_use('ymtrpsnw') ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport 878 IF( iom_use('xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport 879 IF( iom_use('yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport 880 881 DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & 882 & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & 883 & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 884 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 835 836 END DO 837 END DO 838 839 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 840 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 841 & zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. ) 842 843 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) 844 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 845 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 846 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport 847 CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport 848 CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport 849 850 DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 851 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 885 852 886 853 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
