Changeset 11380 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2019-07-31T15:56:02+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/ICE/icedyn_rhg_evp.F90
r11377 r11380 112 112 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 113 113 !! 114 LOGICAL, PARAMETER :: ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 115 ! or only at the main time step (F) 114 116 INTEGER :: ji, jj ! dummy loop indices 115 117 INTEGER :: jter ! local integers … … 135 137 ! 136 138 REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points 137 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV! ice fraction on U/V points139 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 138 140 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 139 141 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 points 143 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 points 140 145 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 stresses 141 147 ! 142 148 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear … … 146 152 ! ! ocean surface (ssh_m) if ice is not embedded 147 153 ! ! ice bottom surface if ice is embedded 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 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 158 159 REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice 159 160 … … 162 163 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 163 164 !! --- diags 164 REAL(wp), DIMENSION(jpi,jpj) :: z msk00165 REAL(wp), DIMENSION(jpi,jpj) :: zswi 165 166 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 166 167 !! --- SIMIP diags 168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice 169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice 170 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 stress 177 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress 167 178 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 168 179 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) … … 253 264 254 265 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 255 IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile256 ELSE ; zkt = 0._wp266 IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN ; zkt = rn_tensile 267 ELSE ; zkt = 0._wp 257 268 ENDIF 258 269 ! … … 297 308 298 309 ! Drag ice-atm. 299 z taux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)300 z tauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)310 zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 311 zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 301 312 302 313 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points … … 305 316 306 317 ! masks 307 zm sk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice308 zm sk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice318 zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 319 zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 309 320 310 321 ! switches 311 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; z msk01x(ji,jj) = 0._wp312 ELSE ; z msk01x(ji,jj) = 1._wp ; ENDIF313 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; z msk01y(ji,jj) = 0._wp314 ELSE ; z msk01y(ji,jj) = 1._wp ; ENDIF322 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zswitchU(ji,jj) = 0._wp 323 ELSE ; zswitchU(ji,jj) = 1._wp ; ENDIF 324 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zswitchV(ji,jj) = 0._wp 325 ELSE ; zswitchV(ji,jj) = 1._wp ; ENDIF 315 326 316 327 END DO … … 328 339 ! ice-bottom stress at U points 329 340 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 330 z taux_base(ji,jj) = -rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )341 zTauU_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 331 342 ! ice-bottom stress at V points 332 343 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 333 z tauy_base(ji,jj) = -rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )344 zTauV_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 334 345 ! ice_bottom stress at T points 335 346 zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 336 tau_icebfr(ji,jj) = -rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )347 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 337 348 END DO 338 349 END DO 339 350 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 340 351 ! 341 ELSE !-- no landfast352 ELSEIF( ln_landfast_home ) THEN !-- Home made 342 353 DO jj = 2, jpjm1 343 354 DO ji = fs_2, fs_jpim1 344 ztaux_base(ji,jj) = 0._wp 345 ztauy_base(ji,jj) = 0._wp 355 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 356 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 357 END DO 358 END DO 359 ! 360 ELSE !-- no landfast 361 DO jj = 2, jpjm1 362 DO ji = fs_2, fs_jpim1 363 zTauU_ib(ji,jj) = 0._wp 364 zTauV_ib(ji,jj) = 0._wp 346 365 END DO 347 366 END DO 348 367 ENDIF 368 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 349 369 350 370 !------------------------------------------------------------------------------! … … 484 504 ! !--- tau_bottom/v_ice 485 505 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 486 zTauB = ztauy_base(ji,jj) / zvel 487 ! !--- OceanBottom-to-Ice stress 488 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 506 zTauB = - zTauV_ib(ji,jj) / zvel 489 507 ! 490 508 ! !--- Coriolis at V-points (energy conserving formulation) 491 zCor V(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &509 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 492 510 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 493 511 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 494 512 ! 495 513 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 496 zTauE = zfV(ji,jj) + z tauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)514 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 497 515 ! 498 516 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 499 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )517 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 500 518 ! 501 519 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 502 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity503 & + zTauE + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)504 &/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast505 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0506 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin507 & ) * zmsk00y(ji,jj)520 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 521 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 522 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 523 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 524 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 525 & ) * zmaskV(ji,jj) 508 526 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 509 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity510 & + zTauE + zTauO * v_ice(ji,jj) ) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)511 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast512 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0513 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin514 & ) * zmsk00y(ji,jj)527 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 528 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 529 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 530 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 531 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 532 & ) * zmaskV(ji,jj) 515 533 ENDIF 516 534 END DO … … 522 540 CALL agrif_interp_ice( 'V' ) 523 541 #endif 524 IF( ln_bdy )CALL bdy_ice_dyn( 'V' )542 IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 525 543 ! 526 544 DO jj = 2, jpjm1 … … 534 552 ! !--- tau_bottom/u_ice 535 553 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 536 zTauB = ztaux_base(ji,jj) / zvel 537 ! !--- OceanBottom-to-Ice stress 538 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 554 zTauB = - zTauU_ib(ji,jj) / zvel 539 555 ! 540 556 ! !--- Coriolis at U-points (energy conserving formulation) 541 zCor U(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &557 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 542 558 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 543 559 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 544 560 ! 545 561 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 546 zTauE = zfU(ji,jj) + z taux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)562 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 547 563 ! 548 564 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 549 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )565 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 550 566 ! 551 567 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 552 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity553 & + zTauE + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)554 &/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast555 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0556 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin557 & ) * zmsk00x(ji,jj)568 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 569 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 570 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 571 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 572 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 573 & ) * zmaskU(ji,jj) 558 574 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 559 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity560 & + zTauE + zTauO * u_ice(ji,jj) ) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)561 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast562 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0563 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin564 & ) * zmsk00x(ji,jj)575 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 576 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 577 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 578 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 579 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 580 & ) * zmaskU(ji,jj) 565 581 ENDIF 566 582 END DO … … 572 588 CALL agrif_interp_ice( 'U' ) 573 589 #endif 574 IF( ln_bdy )CALL bdy_ice_dyn( 'U' )590 IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 575 591 ! 576 592 ELSE ! odd iterations … … 586 602 ! !--- tau_bottom/u_ice 587 603 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 588 zTauB = ztaux_base(ji,jj) / zvel 589 ! !--- OceanBottom-to-Ice stress 590 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 604 zTauB = - zTauU_ib(ji,jj) / zvel 591 605 ! 592 606 ! !--- Coriolis at U-points (energy conserving formulation) 593 zCor U(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &607 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 594 608 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 595 609 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 596 610 ! 597 611 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 598 zTauE = zfU(ji,jj) + z taux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)612 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 599 613 ! 600 614 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 601 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )615 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 602 616 ! 603 617 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 604 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity605 & + zTauE + zTauO * u_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)606 &/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast607 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0608 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin609 & ) * zmsk00x(ji,jj)618 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 619 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 620 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 621 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 622 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 623 & ) * zmaskU(ji,jj) 610 624 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 611 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity612 & + zTauE + zTauO * u_ice(ji,jj) ) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)613 & / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast614 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0615 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin616 & ) * zmsk00x(ji,jj)625 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 626 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 627 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 628 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 629 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 630 & ) * zmaskU(ji,jj) 617 631 ENDIF 618 632 END DO … … 624 638 CALL agrif_interp_ice( 'U' ) 625 639 #endif 626 IF( ln_bdy )CALL bdy_ice_dyn( 'U' )640 IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 627 641 ! 628 642 DO jj = 2, jpjm1 … … 636 650 ! !--- tau_bottom/v_ice 637 651 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 638 zTauB = ztauy_base(ji,jj) / zvel 639 ! !--- OceanBottom-to-Ice stress 640 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 652 zTauB = - zTauV_ib(ji,jj) / zvel 641 653 ! 642 654 ! !--- Coriolis at v-points (energy conserving formulation) 643 zCor V(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &655 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 644 656 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 645 657 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 646 658 ! 647 659 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 648 zTauE = zfV(ji,jj) + z tauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)660 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 649 661 ! 650 662 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 651 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )663 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 652 664 ! 653 665 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 654 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity655 & + zTauE + zTauO * v_ice(ji,jj) )& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)656 &/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast657 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0658 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin659 & ) * zmsk00y(ji,jj)666 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 667 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 668 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 669 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 670 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 671 & ) * zmaskV(ji,jj) 660 672 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 661 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity662 & + zTauE + zTauO * v_ice(ji,jj) ) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)663 & / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast664 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0665 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin666 & ) * zmsk00y(ji,jj)673 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 674 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 675 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 676 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 677 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 678 & ) * zmaskV(ji,jj) 667 679 ENDIF 668 680 END DO … … 674 686 CALL agrif_interp_ice( 'V' ) 675 687 #endif 676 IF( ln_bdy )CALL bdy_ice_dyn( 'V' )688 IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 677 689 ! 678 690 ENDIF … … 689 701 END DO ! end loop over jter ! 690 702 ! ! ==================== ! 703 ! 704 IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN 705 CALL bdy_ice_dyn( 'U' ) 706 CALL bdy_ice_dyn( 'V' ) 707 ENDIF 691 708 ! 692 709 !------------------------------------------------------------------------------! … … 747 764 DO jj = 1, jpj 748 765 DO ji = 1, jpi 749 z msk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice766 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 750 767 END DO 751 768 END DO 752 769 753 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !754 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &755 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN756 !757 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &758 & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )759 !760 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )761 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )762 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )763 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )764 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )765 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )766 ENDIF767 768 770 ! --- divergence, shear and strength --- ! 769 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00) ! divergence770 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00) ! shear771 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) !strength772 773 ! --- stress tensor--- !774 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr')) THEN771 IF( iom_use('icediv') ) CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) ) ! divergence 772 IF( iom_use('iceshe') ) CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) ) ! shear 773 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Ice strength 774 775 ! --- charge ellipse --- ! 776 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 775 777 ! 776 778 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) … … 778 780 DO jj = 2, jpjm1 779 781 DO ji = 2, jpim1 780 zdum1 = ( z msk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point781 & z msk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) &782 & / MAX( 1._wp, z msk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )782 zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 783 & zswi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 784 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 783 785 784 786 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 785 787 786 zdum2 = z msk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )788 zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 787 789 788 790 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) … … 797 799 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 798 800 ! 799 CALL iom_put( 'isig1' , zsig1 ) 800 CALL iom_put( 'isig2' , zsig2 ) 801 CALL iom_put( 'isig3' , zsig3 ) 802 ! 803 ! Stress tensor invariants (normal and shear stress N/m) 804 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 805 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 806 801 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) 802 IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) 803 IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) 804 ! 807 805 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 808 806 ENDIF 809 807 810 808 ! --- SIMIP --- ! 811 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 812 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 813 ! 814 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 815 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 816 817 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 818 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 819 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 820 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 821 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 822 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 823 ENDIF 824 825 IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 826 & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 827 ! 828 ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 829 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 830 ! 809 IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & 810 & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & 811 & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 812 & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN 813 814 ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & 815 & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & 816 & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 817 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) 818 831 819 DO jj = 2, jpjm1 832 820 DO ji = 2, jpim1 821 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 822 823 ! Stress tensor invariants (normal and shear stress N/m) 824 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress 825 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress 826 827 ! Stress terms of the momentum equation (N/m2) 828 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term 829 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 830 831 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term 832 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 833 834 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term 835 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch 836 837 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress 838 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 839 833 840 ! 2D ice mass, snow mass, area transport arrays (X, Y) 834 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)835 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)836 841 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 842 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 843 837 844 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 838 845 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 839 846 840 847 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 841 848 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 842 849 843 850 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 844 851 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 845 846 END DO 847 END DO 848 849 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 850 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 851 & zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. ) 852 853 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) 854 CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport 855 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 856 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport 857 CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport 858 CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport 859 860 DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 861 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 852 853 END DO 854 END DO 855 856 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 857 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 858 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 859 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 860 861 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi , 'U', -1., zdiag_vtau_oi , 'V', -1., & 862 & zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 863 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 864 & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) 865 866 IF( iom_use('normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress 867 IF( iom_use('sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress 868 IF( iom_use('dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) 869 IF( iom_use('dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) 870 IF( iom_use('corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) 871 IF( iom_use('corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) 872 IF( iom_use('intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) 873 IF( iom_use('intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) 874 IF( iom_use('utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) 875 IF( iom_use('vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) 876 IF( iom_use('xmtrpice') ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) 877 IF( iom_use('ymtrpice') ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport 878 IF( iom_use('xmtrpsnw') ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) 879 IF( iom_use('ymtrpsnw') ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport 880 IF( iom_use('xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport 881 IF( iom_use('yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport 882 883 DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & 884 & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & 885 & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 886 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 862 887 863 888 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
