Changeset 11377
- Timestamp:
- 2019-07-31T14:26:56+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization
- Files:
-
- 6 edited
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (2 diffs)
-
cfgs/SHARED/namelist_ice_ref (modified) (1 diff)
-
src/ICE/ice.F90 (modified) (1 diff)
-
src/ICE/icedyn.F90 (modified) (3 diffs)
-
src/ICE/icedyn_rhg.F90 (modified) (1 diff)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/field_def_nemo-ice.xml
r11371 r11377 44 44 <field id="icefrb" long_name="Sea-ice freeboard" standard_name="sea_ice_freeboard" unit="m" /> 45 45 <field id="icealb" long_name="Sea-ice or snow albedo" standard_name="sea_ice_albedo" unit="" detect_missing_value="true" /> 46 <field id="tau_icebfr" long_name="ice friction on ocean bottom for landfast ice" unit="N/2" />47 46 48 47 <!-- melt ponds --> … … 71 70 <field id="utau_oi" long_name="X-component of ocean stress on sea ice" standard_name="sea_ice_base_upward_x_stress" unit="N/m2" /> 72 71 <field id="vtau_oi" long_name="Y-component of ocean stress on sea ice" standard_name="sea_ice_base_upward_y_stress" unit="N/m2" /> 72 <field id="utau_bi" long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress" unit="N/m2" /> 73 <field id="vtau_bi" long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress" unit="N/m2" /> 73 74 <field id="isig1" long_name="1st principal stress component for EVP rhg" unit="" /> 74 75 <field id="isig2" long_name="2nd principal stress component for EVP rhg" unit="" /> -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref
r10911 r11377 56 56 rn_ishlat = 2. ! lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 57 57 ln_landfast_L16 = .false. ! landfast: parameterization from Lemieux 2016 58 ln_landfast_home = .false. ! landfast: parameterization from "home made"59 58 rn_depfra = 0.125 ! fraction of ocean depth that ice must reach to initiate landfast 60 59 ! recommended range: [0.1 ; 0.25] - L16=0.125 - home=0.15 -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90
r11371 r11377 136 136 REAL(wp), PUBLIC :: rn_ishlat !: lateral boundary condition for sea-ice 137 137 LOGICAL , PUBLIC :: ln_landfast_L16 !: landfast ice parameterizationfrom lemieux2016 138 LOGICAL , PUBLIC :: ln_landfast_home !: landfast ice parameterizationfrom home made139 138 REAL(wp), PUBLIC :: rn_depfra !: fraction of ocean depth that ice must reach to initiate landfast ice 140 139 REAL(wp), PUBLIC :: rn_icebfr !: maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn.F90
r11371 r11377 221 221 NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice, & 222 222 & rn_ishlat , & 223 & ln_landfast_L16, ln_landfast_home,rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile223 & ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 224 224 !!------------------------------------------------------------------- 225 225 ! … … 244 244 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 245 245 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 246 WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home247 246 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 248 247 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr … … 271 270 ENDIF 272 271 ! !--- Landfast ice 273 IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home ) tau_icebfr(:,:) = 0._wp 274 ! 275 IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 276 CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 277 ENDIF 272 IF( .NOT.ln_landfast_L16 ) tau_icebfr(:,:) = 0._wp 278 273 ! 279 274 CALL ice_dyn_rdgrft_init ! set ice ridging/rafting parameters -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90
r11317 r11377 69 69 WRITE(numout,*)'ice_dyn_rhg: sea-ice rheology' 70 70 WRITE(numout,*)'~~~~~~~~~~~' 71 ENDIF72 !73 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization74 tau_icebfr(:,:) = 0._wp75 DO jl = 1, jpl76 WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr77 END DO78 71 ENDIF 79 72 ! -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg_evp.F90
r11371 r11377 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 … … 254 253 255 254 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 256 IF( ln_landfast_L16 .OR. ln_landfast_home) THEN ; zkt = rn_tensile257 ELSE 255 IF( ln_landfast_L16 ) THEN ; zkt = rn_tensile 256 ELSE ; zkt = 0._wp 258 257 ENDIF 259 258 ! … … 298 297 299 298 ! Drag ice-atm. 300 z TauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)301 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) 302 301 303 302 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points … … 306 305 307 306 ! masks 308 zm askU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice309 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 310 309 311 310 ! switches 312 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; z switchU(ji,jj) = 0._wp313 ELSE ; z switchU(ji,jj) = 1._wp ; ENDIF314 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; z switchV(ji,jj) = 0._wp315 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 316 315 317 316 END DO … … 329 328 ! ice-bottom stress at U points 330 329 zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 331 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) ) ) 332 331 ! ice-bottom stress at V points 333 332 zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 334 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) ) ) 335 334 ! ice_bottom stress at T points 336 335 zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 337 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) ) ) 338 337 END DO 339 338 END DO 340 339 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 341 340 ! 342 ELSE IF( ln_landfast_home ) THEN !-- Home made341 ELSE !-- no landfast 343 342 DO jj = 2, jpjm1 344 343 DO ji = fs_2, fs_jpim1 345 zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 346 zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 347 END DO 348 END DO 349 ! 350 ELSE !-- no landfast 351 DO jj = 2, jpjm1 352 DO ji = fs_2, fs_jpim1 353 zTauU_ib(ji,jj) = 0._wp 354 zTauV_ib(ji,jj) = 0._wp 344 ztaux_base(ji,jj) = 0._wp 345 ztauy_base(ji,jj) = 0._wp 355 346 END DO 356 347 END DO … … 493 484 ! !--- tau_bottom/v_ice 494 485 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 495 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) 496 489 ! 497 490 ! !--- Coriolis at V-points (energy conserving formulation) 498 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &491 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 499 492 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 500 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) ) ) 501 494 ! 502 495 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 503 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) 504 497 ! 505 498 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 506 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 ) ) ) 507 500 ! 508 501 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 509 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin514 & ) * 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) 515 508 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 516 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity517 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)518 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast519 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0520 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin521 & ) * 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) 522 515 ENDIF 523 516 END DO … … 529 522 CALL agrif_interp_ice( 'V' ) 530 523 #endif 531 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )524 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 532 525 ! 533 526 DO jj = 2, jpjm1 … … 541 534 ! !--- tau_bottom/u_ice 542 535 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 543 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) 544 539 ! 545 540 ! !--- Coriolis at U-points (energy conserving formulation) 546 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &541 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 547 542 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 548 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) ) ) 549 544 ! 550 545 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 551 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) 552 547 ! 553 548 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 554 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 ) ) ) 555 550 ! 556 551 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 557 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity558 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)559 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast560 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0561 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin562 & ) * 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) 563 558 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 564 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity565 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)566 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast567 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0568 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin569 & ) * 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) 570 565 ENDIF 571 566 END DO … … 577 572 CALL agrif_interp_ice( 'U' ) 578 573 #endif 579 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )574 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 580 575 ! 581 576 ELSE ! odd iterations … … 591 586 ! !--- tau_bottom/u_ice 592 587 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 593 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) 594 591 ! 595 592 ! !--- Coriolis at U-points (energy conserving formulation) 596 zCor x(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &593 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 597 594 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 598 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) ) ) 599 596 ! 600 597 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 601 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) 602 599 ! 603 600 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 604 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 ) ) ) 605 602 ! 606 603 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 607 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )& ! previous velocity608 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)609 & )/ MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast610 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )& ! static friction => slow decrease to v=0611 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin612 & ) * 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) 613 610 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 614 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) &! previous velocity615 & + zTauE + zTauO * u_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)616 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast617 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0618 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin619 & ) * 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) 620 617 ENDIF 621 618 END DO … … 627 624 CALL agrif_interp_ice( 'U' ) 628 625 #endif 629 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'U' )626 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 630 627 ! 631 628 DO jj = 2, jpjm1 … … 639 636 ! !--- tau_bottom/v_ice 640 637 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 641 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) 642 641 ! 643 642 ! !--- Coriolis at v-points (energy conserving formulation) 644 zCor y(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &643 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 645 644 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 646 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) ) ) 647 646 ! 648 647 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 649 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) 650 649 ! 651 650 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 652 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 ) ) ) 653 652 ! 654 653 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 655 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )& ! previous velocity656 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)657 & )/ MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast658 &+ ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0659 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin660 & ) * 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) 661 660 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 662 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) &! previous velocity663 & + zTauE + zTauO * v_ice(ji,jj) &! F + tau_ia + Coriolis + spg + tau_io(only ocean part)664 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) &! m/dt + tau_io(only ice part) + landfast665 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) &! static friction => slow decrease to v=0666 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) &! v_ice = v_oce/100 if mass < zmmin & conc < zamin667 & ) * 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) 668 667 ENDIF 669 668 END DO … … 675 674 CALL agrif_interp_ice( 'V' ) 676 675 #endif 677 IF( ln_bdy .AND. ll_bdy_substep )CALL bdy_ice_dyn( 'V' )676 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 678 677 ! 679 678 ENDIF … … 690 689 END DO ! end loop over jter ! 691 690 ! ! ==================== ! 692 !693 IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN694 CALL bdy_ice_dyn( 'U' )695 CALL bdy_ice_dyn( 'V' )696 ENDIF697 691 ! 698 692 !------------------------------------------------------------------------------! … … 757 751 END DO 758 752 759 ! --- ice-ocean and ice-atm. stresses --- ! 760 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') ) THEN 761 CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., & 762 & zTauU_ia, 'U', -1., zTauV_ia, 'V', -1. ) 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 ! 763 760 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 764 761 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 765 CALL iom_put( 'utau_ai' , zTauU_ia * zmsk00 ) 766 CALL iom_put( 'vtau_ai' , zTauV_ia * 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 ) 767 766 ENDIF 768 769 ! --- ice-oceanbottom stress in case of landfast --- ! 770 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr * zmsk00 ) 771 767 772 768 ! --- divergence, shear and strength --- ! 773 769 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence … … 817 813 ! 818 814 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 819 & zCor x, 'U', -1., zCory, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )815 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 820 816 821 817 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 822 818 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 823 CALL iom_put( 'corstrx' , zCor x* zmsk00 ) ! Coriolis force term in force balance (x)824 CALL iom_put( 'corstry' , zCor y* zmsk00 ) ! Coriolis force 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) 825 821 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 826 822 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
Note: See TracChangeset
for help on using the changeset viewer.
