[3614] | 1 | MODULE icbdyn |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE icbdyn *** |
---|
| 4 | !! Iceberg: time stepping routine for iceberg tracking |
---|
| 5 | !!====================================================================== |
---|
[9190] | 6 | !! History : 3.3 ! 2010-01 (Martin&Adcroft) Original code |
---|
| 7 | !! - ! 2011-03 (Madec) Part conversion to NEMO form |
---|
| 8 | !! - ! Removal of mapping from another grid |
---|
| 9 | !! - ! 2011-04 (Alderson) Split into separate modules |
---|
| 10 | !! - ! 2011-05 (Alderson) Replace broken grounding routine with one of |
---|
| 11 | !! - ! Gurvan's suggestions (just like the broken one) |
---|
[3614] | 12 | !!---------------------------------------------------------------------- |
---|
| 13 | USE par_oce ! NEMO parameters |
---|
| 14 | USE dom_oce ! NEMO ocean domain |
---|
| 15 | USE phycst ! NEMO physical constants |
---|
| 16 | ! |
---|
| 17 | USE icb_oce ! define iceberg arrays |
---|
| 18 | USE icbutl ! iceberg utility routines |
---|
| 19 | USE icbdia ! iceberg budget routines |
---|
| 20 | |
---|
| 21 | IMPLICIT NONE |
---|
| 22 | PRIVATE |
---|
| 23 | |
---|
| 24 | PUBLIC icb_dyn ! routine called in icbstp.F90 module |
---|
| 25 | |
---|
| 26 | !!---------------------------------------------------------------------- |
---|
[9598] | 27 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
---|
[5215] | 28 | !! $Id$ |
---|
[10068] | 29 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
[3614] | 30 | !!---------------------------------------------------------------------- |
---|
| 31 | CONTAINS |
---|
| 32 | |
---|
[4990] | 33 | SUBROUTINE icb_dyn( kt ) |
---|
[3614] | 34 | !!---------------------------------------------------------------------- |
---|
| 35 | !! *** ROUTINE icb_dyn *** |
---|
| 36 | !! |
---|
| 37 | !! ** Purpose : iceberg evolution. |
---|
| 38 | !! |
---|
| 39 | !! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010 |
---|
| 40 | !!---------------------------------------------------------------------- |
---|
[9190] | 41 | INTEGER, INTENT(in) :: kt ! |
---|
| 42 | ! |
---|
| 43 | LOGICAL :: ll_bounced |
---|
| 44 | REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 |
---|
| 45 | REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 |
---|
| 46 | REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3 |
---|
| 47 | REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 |
---|
| 48 | REAL(wp) :: zuvel_n, zvvel_n, zxi_n , zyj_n |
---|
| 49 | REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 |
---|
| 50 | TYPE(iceberg), POINTER :: berg |
---|
| 51 | TYPE(point) , POINTER :: pt |
---|
[3614] | 52 | !!---------------------------------------------------------------------- |
---|
[9190] | 53 | ! |
---|
[3614] | 54 | ! 4th order Runge-Kutta to solve: d/dt X = V, d/dt V = A |
---|
| 55 | ! with I.C.'s: X=X1 and V=V1 |
---|
| 56 | ! |
---|
| 57 | ! ; A1=A(X1,V1) |
---|
| 58 | ! X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2) |
---|
| 59 | ! X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3) |
---|
| 60 | ! X4 = X1+ dt*V3 ; V4 = V1+ dt*A3 ; A4=A(X4,V4) |
---|
| 61 | ! |
---|
| 62 | ! Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 |
---|
| 63 | ! Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 |
---|
| 64 | |
---|
| 65 | ! time steps |
---|
| 66 | zdt = berg_dt |
---|
| 67 | zdt_2 = zdt * 0.5_wp |
---|
| 68 | zdt_6 = zdt / 6._wp |
---|
| 69 | |
---|
| 70 | berg => first_berg ! start from the first berg |
---|
| 71 | ! |
---|
| 72 | DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! |
---|
| 73 | ! |
---|
| 74 | pt => berg%current_point |
---|
| 75 | |
---|
[9190] | 76 | ll_bounced = .FALSE. |
---|
[3614] | 77 | |
---|
| 78 | |
---|
| 79 | ! STEP 1 ! |
---|
| 80 | ! ====== ! |
---|
| 81 | zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s |
---|
| 82 | zyj1 = pt%yj ; zvvel1 = pt%vvel |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | ! !** A1 = A(X1,V1) |
---|
| 86 | CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & |
---|
| 87 | & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) |
---|
| 88 | ! |
---|
| 89 | zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt |
---|
| 90 | zv1 = zvvel1 / ze2 |
---|
| 91 | |
---|
| 92 | ! STEP 2 ! |
---|
| 93 | ! ====== ! |
---|
| 94 | ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 |
---|
| 95 | ! position using di/dt & djdt ! V2 in m/s |
---|
| 96 | zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 |
---|
| 97 | zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 |
---|
| 98 | ! |
---|
| 99 | CALL icb_ground( zxi2, zxi1, zu1, & |
---|
[9190] | 100 | & zyj2, zyj1, zv1, ll_bounced ) |
---|
[3614] | 101 | |
---|
| 102 | ! !** A2 = A(X2,V2) |
---|
| 103 | CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & |
---|
| 104 | & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) |
---|
| 105 | ! |
---|
| 106 | zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt |
---|
| 107 | zv2 = zvvel2 / ze2 |
---|
| 108 | ! |
---|
| 109 | ! STEP 3 ! |
---|
| 110 | ! ====== ! |
---|
| 111 | ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) |
---|
| 112 | zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 |
---|
| 113 | zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 |
---|
| 114 | ! |
---|
| 115 | CALL icb_ground( zxi3, zxi1, zu3, & |
---|
[9190] | 116 | & zyj3, zyj1, zv3, ll_bounced ) |
---|
[3614] | 117 | |
---|
| 118 | ! !** A3 = A(X3,V3) |
---|
| 119 | CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & |
---|
| 120 | & zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) |
---|
| 121 | ! |
---|
| 122 | zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt |
---|
| 123 | zv3 = zvvel3 / ze2 |
---|
| 124 | |
---|
| 125 | ! STEP 4 ! |
---|
| 126 | ! ====== ! |
---|
| 127 | ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 |
---|
| 128 | zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 |
---|
| 129 | zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 |
---|
| 130 | |
---|
| 131 | CALL icb_ground( zxi4, zxi1, zu4, & |
---|
[9190] | 132 | & zyj4, zyj1, zv4, ll_bounced ) |
---|
[3614] | 133 | |
---|
| 134 | ! !** A4 = A(X4,V4) |
---|
| 135 | CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & |
---|
| 136 | & zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) |
---|
| 137 | |
---|
| 138 | zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt |
---|
| 139 | zv4 = zvvel4 / ze2 |
---|
| 140 | |
---|
| 141 | ! FINAL STEP ! |
---|
| 142 | ! ========== ! |
---|
| 143 | ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 |
---|
| 144 | ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 |
---|
| 145 | zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 ) |
---|
| 146 | zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 ) |
---|
| 147 | zuvel_n = pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 ) |
---|
| 148 | zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) |
---|
| 149 | |
---|
| 150 | CALL icb_ground( zxi_n, zxi1, zuvel_n, & |
---|
[9190] | 151 | & zyj_n, zyj1, zvvel_n, ll_bounced ) |
---|
[3614] | 152 | |
---|
| 153 | pt%uvel = zuvel_n !** save in berg structure |
---|
| 154 | pt%vvel = zvvel_n |
---|
| 155 | pt%xi = zxi_n |
---|
| 156 | pt%yj = zyj_n |
---|
| 157 | |
---|
| 158 | ! update actual position |
---|
| 159 | pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) |
---|
| 160 | pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) |
---|
| 161 | |
---|
| 162 | berg => berg%next ! switch to the next berg |
---|
| 163 | ! |
---|
| 164 | END DO !== end loop over all bergs ==! |
---|
| 165 | ! |
---|
| 166 | END SUBROUTINE icb_dyn |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | SUBROUTINE icb_ground( pi, pi0, pu, & |
---|
[9190] | 170 | & pj, pj0, pv, ld_bounced ) |
---|
[3614] | 171 | !!---------------------------------------------------------------------- |
---|
| 172 | !! *** ROUTINE icb_ground *** |
---|
| 173 | !! |
---|
| 174 | !! ** Purpose : iceberg grounding. |
---|
| 175 | !! |
---|
| 176 | !! ** Method : - adjust velocity and then put iceberg back to start position |
---|
| 177 | !! NB two possibilities available one of which is hard-coded here |
---|
| 178 | !!---------------------------------------------------------------------- |
---|
| 179 | REAL(wp), INTENT(inout) :: pi , pj ! current iceberg position |
---|
| 180 | REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position |
---|
| 181 | REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities |
---|
| 182 | LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator |
---|
| 183 | ! |
---|
| 184 | INTEGER :: ii, ii0 |
---|
| 185 | INTEGER :: ij, ij0 |
---|
| 186 | INTEGER :: ibounce_method |
---|
| 187 | !!---------------------------------------------------------------------- |
---|
| 188 | ! |
---|
| 189 | ld_bounced = .FALSE. |
---|
| 190 | ! |
---|
| 191 | ii0 = INT( pi0+0.5 ) ; ij0 = INT( pj0+0.5 ) ! initial gridpoint position (T-cell) |
---|
| 192 | ii = INT( pi +0.5 ) ; ij = INT( pj +0.5 ) ! current - - |
---|
| 193 | ! |
---|
| 194 | IF( ii == ii0 .AND. ij == ij0 ) RETURN ! berg remains in the same cell |
---|
| 195 | ! |
---|
| 196 | ! map into current processor |
---|
| 197 | ii0 = mi1( ii0 ) |
---|
| 198 | ij0 = mj1( ij0 ) |
---|
| 199 | ii = mi1( ii ) |
---|
| 200 | ij = mj1( ij ) |
---|
| 201 | ! |
---|
| 202 | IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one |
---|
| 203 | ! |
---|
| 204 | ! From here, berg have reach land: treat grounding/bouncing |
---|
| 205 | ! ------------------------------- |
---|
| 206 | ld_bounced = .TRUE. |
---|
| 207 | |
---|
| 208 | !! not obvious what should happen now |
---|
| 209 | !! if berg tries to enter a land box, the only location we can return it to is the start |
---|
| 210 | !! position (pi0,pj0), since it has to be in a wet box to do any melting; |
---|
| 211 | !! first option is simply to set whole velocity to zero and move back to start point |
---|
| 212 | !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction |
---|
| 213 | !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast |
---|
| 214 | |
---|
| 215 | ibounce_method = 2 |
---|
| 216 | SELECT CASE ( ibounce_method ) |
---|
[9190] | 217 | CASE ( 1 ) |
---|
| 218 | pi = pi0 |
---|
| 219 | pj = pj0 |
---|
| 220 | pu = 0._wp |
---|
| 221 | pv = 0._wp |
---|
| 222 | CASE ( 2 ) |
---|
| 223 | IF( ii0 /= ii ) THEN |
---|
| 224 | pi = pi0 ! return back to the initial position |
---|
| 225 | pu = 0._wp ! zeroing of velocity in the direction of the grounding |
---|
| 226 | ENDIF |
---|
| 227 | IF( ij0 /= ij ) THEN |
---|
| 228 | pj = pj0 ! return back to the initial position |
---|
| 229 | pv = 0._wp ! zeroing of velocity in the direction of the grounding |
---|
| 230 | ENDIF |
---|
[3614] | 231 | END SELECT |
---|
| 232 | ! |
---|
| 233 | END SUBROUTINE icb_ground |
---|
| 234 | |
---|
| 235 | |
---|
| 236 | SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax, & |
---|
| 237 | & pyj, pe2, pvvel, pvvel0, pay, pdt ) |
---|
| 238 | !!---------------------------------------------------------------------- |
---|
| 239 | !! *** ROUTINE icb_accel *** |
---|
| 240 | !! |
---|
| 241 | !! ** Purpose : compute the iceberg acceleration. |
---|
| 242 | !! |
---|
| 243 | !! ** Method : - sum the terms in the momentum budget |
---|
| 244 | !!---------------------------------------------------------------------- |
---|
| 245 | TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg |
---|
| 246 | REAL(wp) , INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential |
---|
| 247 | REAL(wp) , INTENT(in ) :: puvel , pvvel ! berg velocity [m/s] |
---|
| 248 | REAL(wp) , INTENT(in ) :: puvel0, pvvel0 ! initial berg velocity [m/s] |
---|
| 249 | REAL(wp) , INTENT( out) :: pe1, pe2 ! horizontal scale factor at (xi,yj) |
---|
| 250 | REAL(wp) , INTENT(inout) :: pax, pay ! berg acceleration |
---|
| 251 | REAL(wp) , INTENT(in ) :: pdt ! berg time step |
---|
| 252 | ! |
---|
| 253 | REAL(wp), PARAMETER :: pp_alpha = 0._wp ! |
---|
| 254 | REAL(wp), PARAMETER :: pp_beta = 1._wp ! |
---|
| 255 | REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed |
---|
| 256 | REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration |
---|
| 257 | REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! |
---|
| 258 | ! |
---|
| 259 | INTEGER :: itloop |
---|
[9190] | 260 | REAL(wp) :: zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi |
---|
| 261 | REAL(wp) :: zvo, zvi, zva, zvwave, zssh_y |
---|
[3614] | 262 | REAL(wp) :: zff, zT, zD, zW, zL, zM, zF |
---|
| 263 | REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad |
---|
| 264 | REAL(wp) :: z_ocn, z_atm, z_ice |
---|
| 265 | REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop |
---|
| 266 | REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi |
---|
| 267 | REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new |
---|
| 268 | !!---------------------------------------------------------------------- |
---|
| 269 | |
---|
| 270 | ! Interpolate gridded fields to berg |
---|
| 271 | nknberg = berg%number(1) |
---|
| 272 | CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x, & |
---|
| 273 | & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff ) |
---|
| 274 | |
---|
| 275 | zM = berg%current_point%mass |
---|
| 276 | zT = berg%current_point%thickness ! total thickness |
---|
| 277 | zD = ( rn_rho_bergs / pp_rho_seawater ) * zT ! draught (keel depth) |
---|
| 278 | zF = zT - zD ! freeboard |
---|
| 279 | zW = berg%current_point%width |
---|
| 280 | zL = berg%current_point%length |
---|
| 281 | |
---|
| 282 | zhi = MIN( zhi , zD ) |
---|
| 283 | zD_hi = MAX( 0._wp, zD-zhi ) |
---|
| 284 | |
---|
| 285 | ! Wave radiation |
---|
| 286 | zuwave = zua - zuo ; zvwave = zva - zvo ! Use wind speed rel. to ocean for wave model |
---|
| 287 | zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; |
---|
| 288 | ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. |
---|
| 289 | zampl = 0.5 * 0.02025 * zwmod ! This is "a", the wave amplitude |
---|
| 290 | zLwavelength = 0.32 * zwmod ! Surface wave length fitted to data in table at |
---|
| 291 | ! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html |
---|
| 292 | zLcutoff = 0.125 * zLwavelength |
---|
| 293 | zLtop = 0.25 * zLwavelength |
---|
| 294 | zCr = pp_Cr0 * MIN( MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient |
---|
| 295 | ! ! fitted to graph from Carrieres et al., POAC Drift Model. |
---|
| 296 | zwave_rad = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) |
---|
| 297 | zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed |
---|
| 298 | IF( zwmod /= 0._wp ) THEN |
---|
| 299 | zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ? |
---|
| 300 | zvwave = zva/zwmod |
---|
| 301 | ELSE |
---|
| 302 | zuwave = 0. ; zvwave=0. ; zwave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless |
---|
| 303 | ENDIF |
---|
| 304 | |
---|
| 305 | ! Weighted drag coefficients |
---|
| 306 | z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) |
---|
| 307 | z_atm = pp_rho_air / zM * (0.5*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL) |
---|
| 308 | z_ice = pp_rho_ice / zM * (0.5*pp_Cd_iv*zW*zhi ) |
---|
| 309 | IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp |
---|
| 310 | |
---|
| 311 | zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel |
---|
| 312 | ! |
---|
| 313 | DO itloop = 1, 2 ! Iterate on drag coefficients |
---|
| 314 | ! |
---|
| 315 | zus = 0.5 * ( zuveln + puvel ) |
---|
| 316 | zvs = 0.5 * ( zvveln + pvvel ) |
---|
| 317 | zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) |
---|
| 318 | zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) |
---|
| 319 | zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) ) |
---|
| 320 | ! |
---|
| 321 | ! Explicit accelerations |
---|
| 322 | !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & |
---|
| 323 | ! -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) |
---|
| 324 | !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & |
---|
| 325 | ! -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) |
---|
| 326 | zaxe = -grav * zssh_x + zwave_rad * zuwave |
---|
| 327 | zaye = -grav * zssh_y + zwave_rad * zvwave |
---|
| 328 | IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
---|
| 329 | zaxe = zaxe + zff*pvvel0 |
---|
| 330 | zaye = zaye - zff*puvel0 |
---|
| 331 | ELSE |
---|
| 332 | zaxe = zaxe + zff*pvvel |
---|
| 333 | zaye = zaye - zff*puvel |
---|
| 334 | ENDIF |
---|
| 335 | IF( pp_beta > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
---|
| 336 | zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui) |
---|
| 337 | zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi) |
---|
| 338 | ELSE |
---|
| 339 | zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui) |
---|
| 340 | zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi) |
---|
[9190] | 341 | ENDIF |
---|
[3614] | 342 | |
---|
| 343 | ! Solve for implicit accelerations |
---|
| 344 | IF( pp_alpha + pp_beta > 0._wp ) THEN |
---|
| 345 | zlambda = zdrag_ocn + zdrag_atm + zdrag_ice |
---|
| 346 | zA11 = 1._wp + pp_beta *pdt*zlambda |
---|
| 347 | zA12 = pp_alpha*pdt*zff |
---|
| 348 | zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 ) |
---|
| 349 | pax = zdetA * ( zA11*zaxe + zA12*zaye ) |
---|
| 350 | pay = zdetA * ( zA11*zaye - zA12*zaxe ) |
---|
[9190] | 351 | ELSE |
---|
[3614] | 352 | pax = zaxe ; pay = zaye |
---|
[9190] | 353 | ENDIF |
---|
[3614] | 354 | |
---|
| 355 | zuveln = puvel0 + pdt*pax |
---|
| 356 | zvveln = pvvel0 + pdt*pay |
---|
| 357 | ! |
---|
| 358 | END DO ! itloop |
---|
| 359 | |
---|
| 360 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
---|
| 361 | zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg |
---|
| 362 | IF( zspeed > 0._wp ) THEN |
---|
| 363 | zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing |
---|
[9190] | 364 | zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt |
---|
[3614] | 365 | IF( zspeed_new < zspeed ) THEN |
---|
| 366 | zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed |
---|
| 367 | zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction |
---|
| 368 | CALL icb_dia_speed() |
---|
| 369 | ENDIF |
---|
| 370 | ENDIF |
---|
| 371 | ENDIF |
---|
| 372 | ! ! check the speed and acceleration limits |
---|
[10570] | 373 | IF (nn_verbose_level > 0) THEN |
---|
| 374 | IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) & |
---|
| 375 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity' |
---|
| 376 | IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) & |
---|
| 377 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration' |
---|
| 378 | ENDIF |
---|
[3614] | 379 | ! |
---|
| 380 | END SUBROUTINE icb_accel |
---|
| 381 | |
---|
| 382 | !!====================================================================== |
---|
| 383 | END MODULE icbdyn |
---|