Changeset 11413 for NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90
- Timestamp:
- 2019-08-06T17:59:22+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90
r10891 r11413 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)116 114 INTEGER :: ji, jj ! dummy loop indices 117 115 INTEGER :: jter ! local integers … … 137 135 ! 138 136 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 points137 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 140 138 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 141 139 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 140 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 141 ! 148 142 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear … … 152 146 ! ! ocean surface (ssh_m) if ice is not embedded 153 147 ! ! 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 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 159 158 REAL(wp), DIMENSION(jpi,jpj) :: zfmask, zwf ! mask at F points for the ice 160 159 … … 163 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) … … 264 253 265 254 ! 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 255 IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile 256 ELSE ; zkt = 0._wp 268 257 ENDIF 269 258 ! … … 308 297 309 298 ! 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)299 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 300 ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 312 301 313 302 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points … … 316 305 317 306 ! 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 ice307 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 308 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 320 309 321 310 ! switches 322 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; z switchU(ji,jj) = 0._wp323 ELSE ; z switchU(ji,jj) = 1._wp ; ENDIF324 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; z switchV(ji,jj) = 0._wp325 ELSE ; z switchV(ji,jj) = 1._wp ; ENDIF311 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 312 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 313 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 314 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 326 315 327 316 END DO … … 339 328 ! ice-bottom stress at U points 340 329 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 341 z TauU_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )330 ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 342 331 ! ice-bottom stress at V points 343 332 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 344 z TauV_ib(ji,jj) =rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )333 ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 345 334 ! ice_bottom stress at T points 346 335 zvCr = at_i(ji,jj) * rn_depfra * ht_n(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) ) )336 tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 348 337 END DO 349 338 END DO 350 339 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 351 340 ! 352 ELSE IF( ln_landfast_home ) THEN !-- Home made341 ELSE !-- no landfast 353 342 DO jj = 2, jpjm1 354 343 DO ji = fs_2, fs_jpim1 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 344 ztaux_base(ji,jj) = 0._wp 345 ztauy_base(ji,jj) = 0._wp 365 346 END DO 366 347 END DO 367 348 ENDIF 368 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )369 349 370 350 !------------------------------------------------------------------------------! … … 372 352 !------------------------------------------------------------------------------! 373 353 ! 374 ! ! ----------------------!354 ! ! ==================== ! 375 355 DO jter = 1 , nn_nevp ! loop over jter ! 376 ! ! ----------------------!356 ! ! ==================== ! 377 357 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 378 358 ! … … 504 484 ! !--- tau_bottom/v_ice 505 485 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 506 zTauB = - zTauV_ib(ji,jj) / zvel 486 zTauB = ztauy_base(ji,jj) / zvel 487 ! !--- OceanBottom-to-Ice stress 488 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 507 489 ! 508 490 ! !--- Coriolis at V-points (energy conserving formulation) 509 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &491 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 510 492 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 511 493 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 512 494 ! 513 495 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 514 zTauE = zfV(ji,jj) + z TauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)496 zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 515 497 ! 516 498 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 517 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )499 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 518 500 ! 519 501 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 520 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity521 & + 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) + landfast523 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0524 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin525 & ) * zmaskV(ji,jj)502 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 503 & + 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) + landfast 505 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 506 & ) * 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 507 & ) * zmsk00y(ji,jj) 526 508 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 527 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity528 & + 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) + landfast530 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0531 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin532 & ) * zmaskV(ji,jj)509 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 510 & + 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) + landfast 512 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 513 & ) * 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 514 & ) * zmsk00y(ji,jj) 533 515 ENDIF 534 516 END DO … … 540 522 CALL agrif_interp_ice( 'V' ) 541 523 #endif 542 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )524 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 543 525 ! 544 526 DO jj = 2, jpjm1 … … 552 534 ! !--- tau_bottom/u_ice 553 535 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 554 zTauB = - zTauU_ib(ji,jj) / zvel 536 zTauB = ztaux_base(ji,jj) / zvel 537 ! !--- OceanBottom-to-Ice stress 538 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 555 539 ! 556 540 ! !--- Coriolis at U-points (energy conserving formulation) 557 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &541 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 558 542 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 559 543 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 560 544 ! 561 545 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 562 zTauE = zfU(ji,jj) + z TauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)546 zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 563 547 ! 564 548 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 565 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )549 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 566 550 ! 567 551 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 568 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity569 & + 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) + landfast571 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0572 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin573 & ) * zmaskU(ji,jj)552 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 553 & + 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) + landfast 555 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 556 & ) * 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 557 & ) * zmsk00x(ji,jj) 574 558 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 575 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity576 & + 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) + landfast578 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0579 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin580 & ) * zmaskU(ji,jj)559 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 560 & + 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) + landfast 562 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 563 & ) * 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 564 & ) * zmsk00x(ji,jj) 581 565 ENDIF 582 566 END DO … … 588 572 CALL agrif_interp_ice( 'U' ) 589 573 #endif 590 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )574 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 591 575 ! 592 576 ELSE ! odd iterations … … 602 586 ! !--- tau_bottom/u_ice 603 587 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 604 zTauB = - zTauU_ib(ji,jj) / zvel 588 zTauB = ztaux_base(ji,jj) / zvel 589 ! !--- OceanBottom-to-Ice stress 590 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 605 591 ! 606 592 ! !--- Coriolis at U-points (energy conserving formulation) 607 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &593 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 608 594 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 609 595 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 610 596 ! 611 597 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 612 zTauE = zfU(ji,jj) + z TauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)598 zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 613 599 ! 614 600 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 615 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )601 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 616 602 ! 617 603 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 618 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity619 & + 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) + landfast621 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0622 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin623 & ) * zmaskU(ji,jj)604 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 605 & + 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) + landfast 607 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 608 & ) * 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 609 & ) * zmsk00x(ji,jj) 624 610 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 625 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity626 & + 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) + landfast628 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0629 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin630 & ) * zmaskU(ji,jj)611 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 612 & + 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) + landfast 614 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 615 & ) * 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 616 & ) * zmsk00x(ji,jj) 631 617 ENDIF 632 618 END DO … … 638 624 CALL agrif_interp_ice( 'U' ) 639 625 #endif 640 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )626 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 641 627 ! 642 628 DO jj = 2, jpjm1 … … 650 636 ! !--- tau_bottom/v_ice 651 637 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 652 zTauB = - zTauV_ib(ji,jj) / zvel 638 zTauB = ztauy_base(ji,jj) / zvel 639 ! !--- OceanBottom-to-Ice stress 640 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 653 641 ! 654 642 ! !--- Coriolis at v-points (energy conserving formulation) 655 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &643 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 656 644 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 657 645 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 658 646 ! 659 647 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 660 zTauE = zfV(ji,jj) + z TauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)648 zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 661 649 ! 662 650 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 663 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )651 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 664 652 ! 665 653 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 666 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity667 & + 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) + landfast669 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0670 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin671 & ) * zmaskV(ji,jj)654 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 655 & + 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) + landfast 657 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 658 & ) * 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 659 & ) * zmsk00y(ji,jj) 672 660 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 673 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity674 & + 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) + landfast676 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0677 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin678 & ) * zmaskV(ji,jj)661 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 662 & + 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) + landfast 664 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 665 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 666 & ) * zmsk00y(ji,jj) 679 667 ENDIF 680 668 END DO … … 686 674 CALL agrif_interp_ice( 'V' ) 687 675 #endif 688 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )676 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 689 677 ! 690 678 ENDIF … … 701 689 END DO ! end loop over jter ! 702 690 ! ! ==================== ! 703 !704 IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN705 CALL bdy_ice_dyn( 'U' )706 CALL bdy_ice_dyn( 'V' )707 ENDIF708 691 ! 709 692 !------------------------------------------------------------------------------! … … 764 747 DO jj = 1, jpj 765 748 DO ji = 1, jpi 766 z swi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice749 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 767 750 END DO 768 751 END DO 769 752 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') ) THEN 756 ! 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 ENDIF 767 770 768 ! --- divergence, shear and strength --- ! 771 IF( iom_use('icediv') ) CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:)) ! divergence772 IF( iom_use('iceshe') ) CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:)) ! shear773 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Icestrength774 775 ! --- charge ellipse--- !776 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN769 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 770 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 771 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 772 773 ! --- stress tensor --- ! 774 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 777 775 ! 778 776 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) … … 780 778 DO jj = 2, jpjm1 781 779 DO ji = 2, jpim1 782 zdum1 = ( z swi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point783 & z swi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) &784 & / MAX( 1._wp, z swi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )780 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 781 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 782 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 785 783 786 784 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 787 785 788 zdum2 = z swi(ji,jj) / MAX( 1._wp, strength(ji,jj) )786 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 789 787 790 788 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) … … 799 797 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 800 798 ! 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 ! 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 805 807 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 806 808 ENDIF 807 809 808 810 ! --- SIMIP --- ! 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 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 ! 819 831 DO jj = 2, jpjm1 820 832 DO ji = 2, jpim1 821 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice822 823 ! Stress tensor invariants (normal and shear stress N/m)824 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress825 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress826 827 ! Stress terms of the momentum equation (N/m2)828 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term829 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch830 831 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term832 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch833 834 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term835 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch836 837 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress838 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch839 840 833 ! 2D ice mass, snow mass, area transport arrays (X, Y) 841 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch842 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch843 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 844 837 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 845 838 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 846 839 847 840 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 848 841 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 849 842 850 843 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 851 844 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 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 ) 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 ) 887 862 888 863 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
