- Timestamp:
- 2012-04-18T12:42:56+02:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
r3339 r3359 13 13 !! - ! the broken one) 14 14 !!---------------------------------------------------------------------- 15 !!----------------------------------------------------------------------16 !! icb_init : initialise17 !! icb_gen : generate test icebergs18 !! icb_nam : read iceberg namelist19 !!----------------------------------------------------------------------20 15 USE par_oce ! NEMO parameters 21 16 USE dom_oce ! NEMO ocean domain … … 29 24 PRIVATE 30 25 31 PUBLIC evolve_icebergs ! routine called in xxx.F90 module26 PUBLIC evolve_icebergs ! routine called in icbrun.F90 module 32 27 33 28 CONTAINS … … 42 37 !!---------------------------------------------------------------------- 43 38 ! 44 REAL(wp) :: uvel1 , vvel1 , u1, v1, ax1, ay1, xi1 ,yj145 REAL(wp) :: uvel2 , vvel2 , u2, v2, ax2, ay2, xi2 ,yj246 REAL(wp) :: uvel3 , vvel3 , u3, v3, ax3, ay3, xi3 ,yj347 REAL(wp) :: uvel4 , vvel4 , u4, v4, ax4, ay4, xi4 ,yj448 REAL(wp) :: uvel_n, vvel_n , xi_n,yj_n49 REAL(wp) :: zdt, zdt_2, zdt_6, e1,e250 LOGICAL :: bounced39 REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 40 REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 41 REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3 42 REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 43 REAL(wp) :: zuvel_n, zvvel_n , zxi_n, zyj_n 44 REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 45 LOGICAL :: ll_bounced 51 46 TYPE(iceberg), POINTER :: berg 52 47 TYPE(point) , POINTER :: pt … … 77 72 pt => berg%current_point 78 73 79 bounced = .false.74 ll_bounced = .false. 80 75 81 76 82 77 ! STEP 1 ! 83 78 ! ====== ! 84 xi1 = pt%xi ;uvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s85 yj1 = pt%yj ;vvel1 = pt%vvel79 zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s 80 zyj1 = pt%yj ; zvvel1 = pt%vvel 86 81 87 82 88 83 ! !** A1 = A(X1,V1) 89 CALL accel( xi1, e1, uvel1, uvel1,ax1, &90 & berg , yj1, e2, vvel1, vvel1,ay1, zdt_2 )91 ! 92 u1 = uvel1 /e1 !** V1 in d(i,j)/dt93 v1 = vvel1 /e284 CALL accel( zxi1, ze1, zuvel1, zuvel1, zax1, & 85 & berg , zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 86 ! 87 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt 88 zv1 = zvvel1 / ze2 94 89 95 90 ! STEP 2 ! … … 97 92 ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 98 93 ! position using di/dt & djdt ! V2 in m/s 99 xi2 = xi1 + zdt_2 * u1 ; uvel2 = uvel1 + zdt_2 * ax1 100 yj2 = yj1 + zdt_2 * v1 ; vvel2 = vvel1 + zdt_2 * ay1 101 ! 102 CALL adjust_to_ground( xi2, xi1, u1, yj2, yj1, v1, bounced ) 94 zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 95 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 96 ! 97 CALL adjust_to_ground( zxi2, zxi1, zu1, & 98 & zyj2, zyj1, zv1, ll_bounced ) 103 99 104 100 ! !** A2 = A(X2,V2) 105 CALL accel( xi2, e1, uvel2, uvel1,ax2, &106 & berg , yj2, e2, vvel2, vvel1,ay2, zdt_2 )107 ! 108 u2 = uvel2 /e1 !** V2 in d(i,j)/dt109 v2 = vvel2 /e2101 CALL accel( zxi2, ze1, zuvel2, zuvel1, zax2, & 102 & berg , zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 103 ! 104 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt 105 zv2 = zvvel2 / ze2 110 106 ! 111 107 ! STEP 3 ! 112 108 ! ====== ! 113 109 ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) 114 xi3 = xi1 + zdt_2 * u2 ; uvel3 = uvel1 + zdt_2 * ax2 115 yj3 = yj1 + zdt_2 * v2 ; vvel3 = vvel1 + zdt_2 * ay2 116 ! 117 CALL adjust_to_ground( xi3, xi1, u3, yj3, yj1, v3, bounced ) 110 zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 111 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 112 ! 113 CALL adjust_to_ground( zxi3, zxi1, zu3, & 114 & zyj3, zyj1, zv3, ll_bounced ) 118 115 119 116 ! !** A3 = A(X3,V3) 120 CALL accel( xi3, e1, uvel3, uvel1,ax3, &121 & berg , yj3, e2, vvel3, vvel1,ay3, zdt )122 ! 123 u3 = uvel3 /e1 !** V3 in d(i,j)/dt124 v3 = vvel3 /e2117 CALL accel( zxi3, ze1, zuvel3, zuvel1, zax3, & 118 & berg , zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 119 ! 120 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt 121 zv3 = zvvel3 / ze2 125 122 126 123 ! STEP 4 ! 127 124 ! ====== ! 128 125 ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 129 xi4 = xi1 + zdt * u3 ; uvel4 = uvel1 + zdt * ax3 130 yj4 = yj1 + zdt * v3 ; vvel4 = vvel1 + zdt * ay3 131 132 CALL adjust_to_ground( xi4, xi1, u4, yj4, yj1, v4, bounced ) 126 zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 127 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 128 129 CALL adjust_to_ground( zxi4, zxi1, zu4, & 130 & zyj4, zyj1, zv4, ll_bounced ) 133 131 134 132 ! !** A4 = A(X4,V4) 135 CALL accel( xi4, e1, uvel4, uvel1,ax4, &136 & berg , yj4, e2, vvel4, vvel1,ay4, zdt )137 138 u4 = uvel4 /e1 !** V4 in d(i,j)/dt139 v4 = vvel4 /e2133 CALL accel( zxi4, ze1, zuvel4, zuvel1, zax4, & 134 & berg , zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 135 136 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt 137 zv4 = zvvel4 / ze2 140 138 141 139 ! FINAL STEP ! … … 143 141 ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 144 142 ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 145 xi_n = pt%xi + zdt_6 * ( u1 + 2.*(u2 + u3 ) + u4 ) 146 yj_n = pt%yj + zdt_6 * ( v1 + 2.*(v2 + v3 ) + v4 ) 147 uvel_n = pt%uvel + zdt_6 * ( ax1 + 2.*(ax2 + ax3) + ax4 ) 148 vvel_n = pt%vvel + zdt_6 * ( ay1 + 2.*(ay2 + ay3) + ay4 ) 149 150 CALL adjust_to_ground( xi_n, xi1, uvel_n, yj_n, yj1, vvel_n, bounced ) 151 152 pt%uvel = uvel_n !** save in berg structure 153 pt%vvel = vvel_n 154 pt%xi = xi_n 155 pt%yj = yj_n 156 157 ! sga - update actual position 143 zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 ) 144 zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 ) 145 zuvel_n = pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 ) 146 zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) 147 148 CALL adjust_to_ground( zxi_n, zxi1, zuvel_n, & 149 & zyj_n, zyj1, zvvel_n, ll_bounced ) 150 151 pt%uvel = zuvel_n !** save in berg structure 152 pt%vvel = zvvel_n 153 pt%xi = zxi_n 154 pt%yj = zyj_n 155 156 ! update actual position 158 157 pt%lon = bilin_x(glamt, pt%xi, pt%yj ) 159 158 pt%lat = bilin(gphit, pt%xi, pt%yj, 'T', 0, 0 ) … … 166 165 167 166 168 SUBROUTINE adjust_to_ground( pi, pi0, pu, pj, pj0, pv, bounced ) 167 SUBROUTINE adjust_to_ground( pi, pi0, pu, & 168 & pj, pj0, pv, ld_bounced ) 169 169 !!---------------------------------------------------------------------- 170 170 !! *** ROUTINE adjust_to_ground *** … … 177 177 REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position 178 178 REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities 179 LOGICAL , INTENT( out) :: bounced ! bounced indicator 180 ! 181 REAL(wp), PARAMETER :: posn_eps = 0.05_wp ! bouncing distance in (i,j) referential 179 LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator 182 180 ! 183 181 INTEGER :: ii, ii0 184 182 INTEGER :: ij, ij0 185 INTEGER :: bounce_method 186 !!---------------------------------------------------------------------- 187 ! 188 bounced = .false. 189 bounce_method = 2 183 INTEGER :: ibounce_method 184 !!---------------------------------------------------------------------- 185 ! 186 ld_bounced = .FALSE. 190 187 ! 191 188 ii0 = INT( pi0+0.5 ) ; ij0 = INT( pj0+0.5 ) ! initial gridpoint position (T-cell) … … 204 201 ! From here, berg have reach land: treat grounding/bouncing 205 202 ! ------------------------------- 206 bounced = .true.203 ld_bounced = .TRUE. 207 204 208 205 !! not obvious what should happen now … … 213 210 !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast 214 211 215 SELECT CASE ( bounce_method ) 212 ibounce_method = 2 213 SELECT CASE ( ibounce_method ) 216 214 CASE ( 1 ) 217 215 pi = pi0 … … 233 231 234 232 235 SUBROUTINE accel( xi, e1, uvel, uvel0,ax, &236 & berg , yj, e2, vvel, vvel0, ay,dt )233 SUBROUTINE accel( pxi, pe1, puvel, puvel0, pax, & 234 & berg , pyj, pe2, pvvel, pvvel0, pay, pdt ) 237 235 !!---------------------------------------------------------------------- 238 236 !! *** ROUTINE accel *** … … 243 241 !!---------------------------------------------------------------------- 244 242 TYPE(iceberg ), POINTER :: berg 245 REAL(wp), INTENT(in ) :: xi ,yj ! berg position in (i,j) referential246 REAL(wp), INTENT(in ) :: uvel ,vvel ! berg velocity [m/s]247 REAL(wp), INTENT(in ) :: uvel0,vvel0 ! initial berg velocity [m/s]248 REAL(wp), INTENT( out) :: e1,e2 ! horizontal scale factor at (xi,yj)249 REAL(wp), INTENT(inout) :: ax,ay ! berg acceleration250 REAL(wp), INTENT(in ) :: dt ! berg time step251 ! 252 REAL(wp), PARAMETER :: alpha = 0._wp !253 REAL(wp), PARAMETER :: beta = 1._wp !254 REAL(wp), PARAMETER :: vel_lim =15._wp ! max allowed berg speed255 REAL(wp), PARAMETER :: accel_lim = 1.e-2_wp ! max allowed berg acceleration256 REAL(wp), PARAMETER :: Cr0 = 0.06_wp !243 REAL(wp), INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential 244 REAL(wp), INTENT(in ) :: puvel , pvvel ! berg velocity [m/s] 245 REAL(wp), INTENT(in ) :: puvel0, pvvel0 ! initial berg velocity [m/s] 246 REAL(wp), INTENT( out) :: pe1, pe2 ! horizontal scale factor at (xi,yj) 247 REAL(wp), INTENT(inout) :: pax, pay ! berg acceleration 248 REAL(wp), INTENT(in ) :: pdt ! berg time step 249 ! 250 REAL(wp), PARAMETER :: pp_alpha = 0._wp ! 251 REAL(wp), PARAMETER :: pp_beta = 1._wp ! 252 REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed 253 REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration 254 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 257 255 ! 258 256 INTEGER :: itloop 259 REAL(wp) :: uo, vo, ui, vi, ua, va, uwave, vwave, ssh_x, ssh_y, sst, cn,hi260 REAL(wp) :: zff, T, D, W, L, M,F261 REAL(wp) :: drag_ocn, drag_atm, drag_ice,wave_rad262 REAL(wp) :: c_ocn, c_atm, c_ice263 REAL(wp) :: ampl, wmod, Cr, Lwavelength, Lcutoff,Ltop264 REAL(wp) :: lambda, detA, A11, A12, axe, aye,D_hi265 REAL(wp) :: uveln, vveln, us, vs, speed, loc_dx, new_speed257 REAL(wp) :: zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi 258 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 259 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 260 REAL(wp) :: z_ocn, z_atm, z_ice 261 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 262 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 263 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 266 264 !!---------------------------------------------------------------------- 267 265 268 266 ! Interpolate gridded fields to berg 269 knberg = berg%number(1)270 CALL interp_flds( xi, e1, uo, ui, ua,ssh_x, &271 & yj, e2, vo, vi, va, ssh_y, sst, cn,hi, zff )272 273 M = berg%current_point%mass274 T = berg%current_point%thickness ! total thickness275 D = ( rn_rho_bergs / rho_seawater ) *T ! draught (keel depth)276 F = T -D ! freeboard277 W = berg%current_point%width278 L = berg%current_point%length279 280 hi = MIN( hi ,D )281 D_hi = MAX( 0._wp, D-hi )267 nknberg = berg%number(1) 268 CALL interp_flds( pxi, pe1, zuo, zui, zua, zssh_x, & 269 & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff ) 270 271 zM = berg%current_point%mass 272 zT = berg%current_point%thickness ! total thickness 273 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT ! draught (keel depth) 274 zF = zT - zD ! freeboard 275 zW = berg%current_point%width 276 zL = berg%current_point%length 277 278 zhi = MIN( zhi , zD ) 279 zD_hi = MAX( 0._wp, zD-zhi ) 282 280 283 281 ! Wave radiation 284 uwave = ua - uo ; vwave = va -vo ! Use wind speed rel. to ocean for wave model285 wmod = uwave*uwave + vwave*vwave ! The wave amplitude and length depend on the current;282 zuwave = zua - zuo ; zvwave = zva - zvo ! Use wind speed rel. to ocean for wave model 283 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 286 284 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. 287 ampl = 0.5 * 0.02025 *wmod ! This is "a", the wave amplitude288 Lwavelength = 0.32 *wmod ! Surface wave length fitted to data in table at285 zampl = 0.5 * 0.02025 * zwmod ! This is "a", the wave amplitude 286 zLwavelength = 0.32 * zwmod ! Surface wave length fitted to data in table at 289 287 ! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html 290 Lcutoff = 0.125 *Lwavelength291 Ltop = 0.25 *Lwavelength292 Cr = Cr0 * MIN( MAX( 0., (L-Lcutoff) / ((Ltop-Lcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient288 zLcutoff = 0.125 * zLwavelength 289 zLtop = 0.25 * zLwavelength 290 zCr = pp_Cr0 * MIN( MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient 293 291 ! ! fitted to graph from Carrieres et al., POAC Drift Model. 294 wave_rad = 0.5 * rho_seawater / M * Cr * grav * ampl * MIN( ampl,F ) * (2.*W*L) / (W+L)295 wmod = SQRT( ua*ua + va*va ) ! Wind speed296 IF( wmod /= 0._wp ) THEN297 uwave = ua/wmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ?298 vwave = va/wmod292 zwave_rad = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) 293 zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed 294 IF( zwmod /= 0._wp ) THEN 295 zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ? 296 zvwave = zva/zwmod 299 297 ELSE 300 uwave = 0. ; vwave=0. ;wave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless298 zuwave = 0. ; zvwave=0. ; zwave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless 301 299 ENDIF 302 300 303 301 ! Weighted drag coefficients 304 c_ocn = rho_seawater / M * (0.5*Cd_wv*W*(D_hi)+Cd_wh*W*L)305 c_atm = rho_air / M * (0.5*Cd_av*W*F +Cd_ah*W*L)306 c_ice = rho_ice / M * (0.5*Cd_iv*W*hi )307 IF( abs( ui) + abs(vi) == 0._wp ) c_ice = 0._wp308 309 uveln = uvel ; vveln =vvel ! Copy starting uvel, vvel302 z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) 303 z_atm = pp_rho_air / zM * (0.5*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL) 304 z_ice = pp_rho_ice / zM * (0.5*pp_Cd_iv*zW*zhi ) 305 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 306 307 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 310 308 ! 311 309 DO itloop = 1, 2 ! Iterate on drag coefficients 312 310 ! 313 us = 0.5 * ( uveln +uvel )314 vs = 0.5 * ( vveln +vvel )315 drag_ocn = c_ocn * SQRT( (us-uo)*(us-uo) + (vs-vo)*(vs-vo) )316 drag_atm = c_atm * SQRT( (us-ua)*(us-ua) + (vs-va)*(vs-va) )317 drag_ice = c_ice * SQRT( (us-ui)*(us-ui) + (vs-vi)*(vs-vi) )311 zus = 0.5 * ( zuveln + puvel ) 312 zvs = 0.5 * ( zvveln + pvvel ) 313 zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) 314 zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) 315 zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) ) 318 316 ! 319 317 ! Explicit accelerations 320 ! axe= zff*vvel -grav*ssh_x +wave_rad*uwave &321 ! - drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)322 ! aye=-zff*uvel -grav*ssh_y +wave_rad*vwave &323 ! - drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)324 axe = -grav * ssh_x + wave_rad *uwave325 aye = -grav * ssh_y + wave_rad *vwave326 IF( alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest327 axe = axe + zff*vvel0328 aye = aye - zff*uvel0318 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 319 ! -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 320 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 321 ! -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 322 zaxe = -grav * zssh_x + zwave_rad * zuwave 323 zaye = -grav * zssh_y + zwave_rad * zvwave 324 IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest 325 zaxe = zaxe + zff*pvvel0 326 zaye = zaye - zff*puvel0 329 327 else 330 axe=axe+zff*vvel331 aye=aye-zff*uvel328 zaxe=zaxe+zff*pvvel 329 zaye=zaye-zff*puvel 332 330 endif 333 if ( beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest334 axe=axe-drag_ocn*(uvel0-uo) -drag_atm*(uvel0-ua) -drag_ice*(uvel0-ui)335 aye=aye-drag_ocn*(vvel0-vo) -drag_atm*(vvel0-va) -drag_ice*(vvel0-vi)331 if (pp_beta>0.) then ! If implicit, use time-level (n) rather than RK4 latest 332 zaxe=zaxe-zdrag_ocn*(puvel0-zuo) -zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui) 333 zaye=zaye-zdrag_ocn*(pvvel0-zvo) -zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi) 336 334 else 337 axe=axe-drag_ocn*(uvel-uo) -drag_atm*(uvel-ua) -drag_ice*(uvel-ui)338 aye=aye-drag_ocn*(vvel-vo) -drag_atm*(vvel-va) -drag_ice*(vvel-vi)335 zaxe=zaxe-zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 336 zaye=zaye-zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 339 337 endif 340 338 341 339 ! Solve for implicit accelerations 342 IF( alpha +beta > 0._wp ) THEN343 lambda = drag_ocn + drag_atm +drag_ice344 A11 = 1.+beta*dt*lambda345 A12 = alpha*dt*zff346 detA = 1._wp / ( A11*A11 + A12*A12 )347 ax = detA * ( A11*axe+A12*aye )348 ay = detA * ( A11*aye-A12*axe )340 IF( pp_alpha + pp_beta > 0._wp ) THEN 341 zlambda = zdrag_ocn + zdrag_atm + zdrag_ice 342 zA11 = 1.+pp_beta*pdt*zlambda 343 zA12 = pp_alpha*pdt*zff 344 zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 ) 345 pax = zdetA * ( zA11*zaxe+zA12*zaye ) 346 pay = zdetA * ( zA11*zaye-zA12*zaxe ) 349 347 else 350 ax=axe ; ay=aye348 pax=zaxe ; pay=zaye 351 349 endif 352 350 353 uveln = uvel0 + dt*ax354 vveln = vvel0 + dt*ay351 zuveln = puvel0 + pdt*pax 352 zvveln = pvvel0 + pdt*pay 355 353 ! 356 354 END DO ! itloop 357 355 358 356 IF( rn_speed_limit > 0.) THEN ! Limit speed of bergs based on a CFL criteria (if asked) 359 speed = SQRT( uveln*uveln + vveln*vveln ) ! Speed of berg360 IF( speed > 0.) THEN361 loc_dx = MIN( e1,e2 ) ! minimum grid spacing362 new_speed = loc_dx /dt * rn_speed_limit ! Speed limit as a factor of dx / dt363 IF( new_speed <speed ) THEN364 uveln = uveln * ( new_speed /speed ) ! Scale velocity to reduce speed365 vveln = vveln * ( new_speed /speed ) ! without changing the direction357 zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg 358 IF( zspeed > 0.) THEN 359 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 360 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt 361 IF( zspeed_new < zspeed ) THEN 362 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed 363 zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction 366 364 CALL speed_budget() 367 365 ENDIF … … 369 367 ENDIF 370 368 ! ! check the speed and acceleration limits 371 IF( ABS( uveln ) > vel_lim .OR. ABS( vveln ) >vel_lim ) &369 IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) & 372 370 WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity' 373 IF( ABS( ax ) > accel_lim .OR. ABS( ay ) >accel_lim ) &371 IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) & 374 372 WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration' 375 373 !
Note: See TracChangeset
for help on using the changeset viewer.
