| 1 | MODULE icbdyn |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE icbdyn *** |
|---|
| 4 | !! Iceberg: time stepping routine for iceberg tracking |
|---|
| 5 | !!====================================================================== |
|---|
| 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) |
|---|
| 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 | !!---------------------------------------------------------------------- |
|---|
| 27 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
|---|
| 28 | !! $Id$ |
|---|
| 29 | !! Software governed by the CeCILL license (see ./LICENSE) |
|---|
| 30 | !!---------------------------------------------------------------------- |
|---|
| 31 | CONTAINS |
|---|
| 32 | |
|---|
| 33 | SUBROUTINE icb_dyn( kt ) |
|---|
| 34 | !!---------------------------------------------------------------------- |
|---|
| 35 | !! *** ROUTINE icb_dyn *** |
|---|
| 36 | !! |
|---|
| 37 | !! ** Purpose : iceberg evolution. |
|---|
| 38 | !! |
|---|
| 39 | !! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010 |
|---|
| 40 | !!---------------------------------------------------------------------- |
|---|
| 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 |
|---|
| 52 | !!---------------------------------------------------------------------- |
|---|
| 53 | ! |
|---|
| 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 | |
|---|
| 76 | ll_bounced = .FALSE. |
|---|
| 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, 1 ) |
|---|
| 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, & |
|---|
| 100 | & zyj2, zyj1, zv1, ll_bounced ) |
|---|
| 101 | |
|---|
| 102 | ! !** A2 = A(X2,V2) |
|---|
| 103 | CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & |
|---|
| 104 | & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 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, & |
|---|
| 116 | & zyj3, zyj1, zv3, ll_bounced ) |
|---|
| 117 | |
|---|
| 118 | ! !** A3 = A(X3,V3) |
|---|
| 119 | CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & |
|---|
| 120 | & zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 3 ) |
|---|
| 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, & |
|---|
| 132 | & zyj4, zyj1, zv4, ll_bounced ) |
|---|
| 133 | |
|---|
| 134 | ! !** A4 = A(X4,V4) |
|---|
| 135 | CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & |
|---|
| 136 | & zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 4 ) |
|---|
| 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, & |
|---|
| 151 | & zyj_n, zyj1, zvvel_n, ll_bounced ) |
|---|
| 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, & |
|---|
| 170 | & pj, pj0, pv, ld_bounced ) |
|---|
| 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 ) |
|---|
| 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 |
|---|
| 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, pstage ) |
|---|
| 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 | INTEGER , INTENT(in ) :: pstage ! RK4 stage (for speed limiter) |
|---|
| 253 | ! |
|---|
| 254 | REAL(wp), PARAMETER :: pp_alpha = 0._wp ! |
|---|
| 255 | REAL(wp), PARAMETER :: pp_beta = 1._wp ! |
|---|
| 256 | REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed |
|---|
| 257 | REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration |
|---|
| 258 | REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! |
|---|
| 259 | ! |
|---|
| 260 | INTEGER :: itloop |
|---|
| 261 | REAL(wp) :: zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss |
|---|
| 262 | REAL(wp) :: zvo, zvi, zva, zvwave, zssh_y |
|---|
| 263 | REAL(wp) :: zff, zT, zD, zW, zL, zM, zF |
|---|
| 264 | REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad |
|---|
| 265 | REAL(wp) :: z_ocn, z_atm, z_ice |
|---|
| 266 | REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop |
|---|
| 267 | REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi |
|---|
| 268 | REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale |
|---|
| 269 | !!---------------------------------------------------------------------- |
|---|
| 270 | |
|---|
| 271 | ! Interpolate gridded fields to berg |
|---|
| 272 | nknberg = berg%number(1) |
|---|
| 273 | CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x, & |
|---|
| 274 | & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss ) |
|---|
| 275 | |
|---|
| 276 | zM = berg%current_point%mass |
|---|
| 277 | zT = berg%current_point%thickness ! total thickness |
|---|
| 278 | zD = ( rn_rho_bergs / pp_rho_seawater ) * zT ! draught (keel depth) |
|---|
| 279 | zF = zT - zD ! freeboard |
|---|
| 280 | zW = berg%current_point%width |
|---|
| 281 | zL = berg%current_point%length |
|---|
| 282 | |
|---|
| 283 | zhi = MIN( zhi , zD ) |
|---|
| 284 | zD_hi = MAX( 0._wp, zD-zhi ) |
|---|
| 285 | |
|---|
| 286 | ! Wave radiation |
|---|
| 287 | zuwave = zua - zuo ; zvwave = zva - zvo ! Use wind speed rel. to ocean for wave model |
|---|
| 288 | zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; |
|---|
| 289 | ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. |
|---|
| 290 | zampl = 0.5 * 0.02025 * zwmod ! This is "a", the wave amplitude |
|---|
| 291 | zLwavelength = 0.32 * zwmod ! Surface wave length fitted to data in table at |
|---|
| 292 | ! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html |
|---|
| 293 | zLcutoff = 0.125 * zLwavelength |
|---|
| 294 | zLtop = 0.25 * zLwavelength |
|---|
| 295 | zCr = pp_Cr0 * MIN( MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient |
|---|
| 296 | ! ! fitted to graph from Carrieres et al., POAC Drift Model. |
|---|
| 297 | zwave_rad = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) |
|---|
| 298 | zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed |
|---|
| 299 | IF( zwmod /= 0._wp ) THEN |
|---|
| 300 | zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ? |
|---|
| 301 | zvwave = zva/zwmod |
|---|
| 302 | ELSE |
|---|
| 303 | zuwave = 0. ; zvwave=0. ; zwave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless |
|---|
| 304 | ENDIF |
|---|
| 305 | |
|---|
| 306 | ! Weighted drag coefficients |
|---|
| 307 | z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) |
|---|
| 308 | z_atm = pp_rho_air / zM * (0.5*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL) |
|---|
| 309 | z_ice = pp_rho_ice / zM * (0.5*pp_Cd_iv*zW*zhi ) |
|---|
| 310 | IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp |
|---|
| 311 | |
|---|
| 312 | zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel |
|---|
| 313 | ! |
|---|
| 314 | DO itloop = 1, 2 ! Iterate on drag coefficients |
|---|
| 315 | ! |
|---|
| 316 | zus = 0.5 * ( zuveln + puvel ) |
|---|
| 317 | zvs = 0.5 * ( zvveln + pvvel ) |
|---|
| 318 | zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) |
|---|
| 319 | zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) |
|---|
| 320 | zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) ) |
|---|
| 321 | ! |
|---|
| 322 | ! Explicit accelerations |
|---|
| 323 | !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & |
|---|
| 324 | ! -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) |
|---|
| 325 | !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & |
|---|
| 326 | ! -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) |
|---|
| 327 | zaxe = -grav * zssh_x + zwave_rad * zuwave |
|---|
| 328 | zaye = -grav * zssh_y + zwave_rad * zvwave |
|---|
| 329 | IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
|---|
| 330 | zaxe = zaxe + zff*pvvel0 |
|---|
| 331 | zaye = zaye - zff*puvel0 |
|---|
| 332 | ELSE |
|---|
| 333 | zaxe = zaxe + zff*pvvel |
|---|
| 334 | zaye = zaye - zff*puvel |
|---|
| 335 | ENDIF |
|---|
| 336 | IF( pp_beta > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
|---|
| 337 | zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui) |
|---|
| 338 | zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi) |
|---|
| 339 | ELSE |
|---|
| 340 | zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui) |
|---|
| 341 | zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi) |
|---|
| 342 | ENDIF |
|---|
| 343 | |
|---|
| 344 | ! Solve for implicit accelerations |
|---|
| 345 | IF( pp_alpha + pp_beta > 0._wp ) THEN |
|---|
| 346 | zlambda = zdrag_ocn + zdrag_atm + zdrag_ice |
|---|
| 347 | zA11 = 1._wp + pp_beta *pdt*zlambda |
|---|
| 348 | zA12 = pp_alpha*pdt*zff |
|---|
| 349 | zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 ) |
|---|
| 350 | pax = zdetA * ( zA11*zaxe + zA12*zaye ) |
|---|
| 351 | pay = zdetA * ( zA11*zaye - zA12*zaxe ) |
|---|
| 352 | ELSE |
|---|
| 353 | pax = zaxe ; pay = zaye |
|---|
| 354 | ENDIF |
|---|
| 355 | |
|---|
| 356 | zuveln = puvel0 + pdt*pax |
|---|
| 357 | zvveln = pvvel0 + pdt*pay |
|---|
| 358 | ! |
|---|
| 359 | END DO ! itloop |
|---|
| 360 | |
|---|
| 361 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
|---|
| 362 | zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg |
|---|
| 363 | IF( zspeed > 0._wp ) THEN |
|---|
| 364 | zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing |
|---|
| 365 | ! cfl scale is function of the RK4 step |
|---|
| 366 | IF( pstage .le. 2 ) THEN |
|---|
| 367 | zcfl_scale=0.5 |
|---|
| 368 | ELSE |
|---|
| 369 | zcfl_scale=1.0 |
|---|
| 370 | ENDIF |
|---|
| 371 | zspeed_new = zloc_dx / pdt * rn_speed_limit * zcfl_scale ! Speed limit as a factor of dx / dt |
|---|
| 372 | IF( zspeed_new < zspeed ) THEN |
|---|
| 373 | zvel_factor = zspeed_new / zspeed |
|---|
| 374 | zuveln = zuveln * ( zvel_factor ) ! Scale velocity to reduce speed |
|---|
| 375 | zvveln = zvveln * ( zvel_factor ) ! without changing the direction |
|---|
| 376 | pax = (zuveln - puvel0)/pdt |
|---|
| 377 | pay = (zvveln - pvvel0)/pdt |
|---|
| 378 | WRITE(numicb,*) 'speeding ticket (zspeed_new/zspeed): ',zvel_factor, pdt, zcfl_scale |
|---|
| 379 | CALL FLUSH(numicb) |
|---|
| 380 | CALL icb_dia_speed(zvel_factor,pstage) |
|---|
| 381 | ENDIF |
|---|
| 382 | ENDIF |
|---|
| 383 | ENDIF |
|---|
| 384 | ! ! check the speed and acceleration limits |
|---|
| 385 | IF (nn_verbose_level > 0) THEN |
|---|
| 386 | IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) & |
|---|
| 387 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity' |
|---|
| 388 | IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) & |
|---|
| 389 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration' |
|---|
| 390 | ENDIF |
|---|
| 391 | ! |
|---|
| 392 | END SUBROUTINE icb_accel |
|---|
| 393 | |
|---|
| 394 | !!====================================================================== |
|---|
| 395 | END MODULE icbdyn |
|---|
