| 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, l_loop |
|---|
| 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) :: zvel_factor, zuvel_rk4, zvvel_rk4 |
|---|
| 50 | REAL(wp) :: zvel_factor_keep(100) |
|---|
| 51 | REAL(wp) :: zspeed, zspeed_new, zloc_dx |
|---|
| 52 | REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 |
|---|
| 53 | INTEGER :: icount |
|---|
| 54 | TYPE(iceberg), POINTER :: berg |
|---|
| 55 | TYPE(point) , POINTER :: pt |
|---|
| 56 | !!---------------------------------------------------------------------- |
|---|
| 57 | ! |
|---|
| 58 | ! 4th order Runge-Kutta to solve: d/dt X = V, d/dt V = A |
|---|
| 59 | ! with I.C.'s: X=X1 and V=V1 |
|---|
| 60 | ! |
|---|
| 61 | ! ; A1=A(X1,V1) |
|---|
| 62 | ! X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 ; A2=A(X2,V2) |
|---|
| 63 | ! X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2 ; A3=A(X3,V3) |
|---|
| 64 | ! X4 = X1+ dt*V3 ; V4 = V1+ dt*A3 ; A4=A(X4,V4) |
|---|
| 65 | ! |
|---|
| 66 | ! Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 |
|---|
| 67 | ! Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 |
|---|
| 68 | |
|---|
| 69 | ! time steps |
|---|
| 70 | zdt = berg_dt |
|---|
| 71 | zdt_2 = zdt * 0.5_wp |
|---|
| 72 | zdt_6 = zdt / 6._wp |
|---|
| 73 | |
|---|
| 74 | berg => first_berg ! start from the first berg |
|---|
| 75 | ! |
|---|
| 76 | DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! |
|---|
| 77 | ! |
|---|
| 78 | pt => berg%current_point |
|---|
| 79 | |
|---|
| 80 | ll_bounced = .FALSE. |
|---|
| 81 | |
|---|
| 82 | zvel_factor = 1.0 |
|---|
| 83 | zvel_factor_keep(:) = 0.0 |
|---|
| 84 | l_loop = .true. |
|---|
| 85 | icount = 0 |
|---|
| 86 | |
|---|
| 87 | DO WHILE(l_loop) |
|---|
| 88 | |
|---|
| 89 | l_loop = .false. |
|---|
| 90 | icount=icount+1 |
|---|
| 91 | zvel_factor_keep(icount) = zvel_factor |
|---|
| 92 | |
|---|
| 93 | ! STEP 1 ! |
|---|
| 94 | ! ====== ! |
|---|
| 95 | zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s |
|---|
| 96 | zyj1 = pt%yj ; zvvel1 = pt%vvel |
|---|
| 97 | ! !** A1 = A(X1,V1) |
|---|
| 98 | CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & |
|---|
| 99 | & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) |
|---|
| 100 | ! |
|---|
| 101 | zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt |
|---|
| 102 | zv1 = zvvel1 / ze2 |
|---|
| 103 | |
|---|
| 104 | ! STEP 2 ! |
|---|
| 105 | ! ====== ! |
|---|
| 106 | ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 |
|---|
| 107 | ! position using di/dt & djdt ! V2 in m/s |
|---|
| 108 | zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zvel_factor * ( zuvel1 + zdt_2 * zax1 ) |
|---|
| 109 | zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvel_factor * ( zvvel1 + zdt_2 * zay1 ) |
|---|
| 110 | ! |
|---|
| 111 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
|---|
| 112 | CALL icb_speed_limit(ze1, ze2, zdt_2, zuvel2, zvvel2, rn_speed_limit, l_loop, zvel_factor, icount, 1) |
|---|
| 113 | IF(l_loop) CYCLE |
|---|
| 114 | ENDIF |
|---|
| 115 | ! |
|---|
| 116 | CALL icb_ground( zxi2, zxi1, zu1, & |
|---|
| 117 | & zyj2, zyj1, zv1, ll_bounced ) |
|---|
| 118 | |
|---|
| 119 | ! !** A2 = A(X2,V2) |
|---|
| 120 | CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & |
|---|
| 121 | & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) |
|---|
| 122 | ! |
|---|
| 123 | zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt |
|---|
| 124 | zv2 = zvvel2 / ze2 |
|---|
| 125 | ! |
|---|
| 126 | ! STEP 3 ! |
|---|
| 127 | ! ====== ! |
|---|
| 128 | ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) |
|---|
| 129 | zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zvel_factor * ( zuvel1 + zdt_2 * zax2 ) |
|---|
| 130 | zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvel_factor * ( zvvel1 + zdt_2 * zay2 ) |
|---|
| 131 | ! |
|---|
| 132 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
|---|
| 133 | CALL icb_speed_limit(ze1, ze2, zdt, zuvel3, zvvel3, rn_speed_limit, l_loop, zvel_factor, icount, 2) |
|---|
| 134 | IF(l_loop) CYCLE |
|---|
| 135 | ENDIF |
|---|
| 136 | ! |
|---|
| 137 | CALL icb_ground( zxi3, zxi1, zu3, & |
|---|
| 138 | & zyj3, zyj1, zv3, ll_bounced ) |
|---|
| 139 | |
|---|
| 140 | ! !** A3 = A(X3,V3) |
|---|
| 141 | CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & |
|---|
| 142 | & zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) |
|---|
| 143 | ! |
|---|
| 144 | zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt |
|---|
| 145 | zv3 = zvvel3 / ze2 |
|---|
| 146 | |
|---|
| 147 | ! STEP 4 ! |
|---|
| 148 | ! ====== ! |
|---|
| 149 | ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 |
|---|
| 150 | zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zvel_factor * ( zuvel1 + zdt * zax3 ) |
|---|
| 151 | zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvel_factor * ( zvvel1 + zdt * zay3 ) |
|---|
| 152 | ! |
|---|
| 153 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
|---|
| 154 | zuvel_rk4 = zuvel1 + 2.*(zuvel2 + zuvel3) + zuvel4 |
|---|
| 155 | zvvel_rk4 = zvvel1 + 2.*(zvvel2 + zvvel3) + zuvel4 |
|---|
| 156 | CALL icb_speed_limit(ze1, ze2, zdt_6, zuvel_rk4, zvvel_rk4, rn_speed_limit, l_loop, zvel_factor, icount, 3) |
|---|
| 157 | IF(l_loop) CYCLE |
|---|
| 158 | ENDIF |
|---|
| 159 | |
|---|
| 160 | CALL icb_ground( zxi4, zxi1, zu4, & |
|---|
| 161 | & zyj4, zyj1, zv4, ll_bounced ) |
|---|
| 162 | |
|---|
| 163 | ! !** A4 = A(X4,V4) |
|---|
| 164 | CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & |
|---|
| 165 | & zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) |
|---|
| 166 | |
|---|
| 167 | zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt |
|---|
| 168 | zv4 = zvvel4 / ze2 |
|---|
| 169 | |
|---|
| 170 | ! FINAL STEP ! |
|---|
| 171 | ! ========== ! |
|---|
| 172 | ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 |
|---|
| 173 | ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 |
|---|
| 174 | zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 ) |
|---|
| 175 | zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 ) |
|---|
| 176 | zuvel_n = zvel_factor * ( pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 ) ) |
|---|
| 177 | zvvel_n = zvel_factor * ( pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) ) |
|---|
| 178 | |
|---|
| 179 | IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked) |
|---|
| 180 | CALL icb_speed_limit(ze1, ze2, zdt_2, zuvel_n, zvvel_n, rn_speed_limit, l_loop, zvel_factor, icount, 4) |
|---|
| 181 | ENDIF |
|---|
| 182 | |
|---|
| 183 | END DO |
|---|
| 184 | |
|---|
| 185 | IF( icount > 5 ) WRITE(numicb, '("Large icount : vel_factors = ",i4,":",30(",",f12.8))') icount,zvel_factor_keep(1:icount) |
|---|
| 186 | |
|---|
| 187 | CALL icb_ground( zxi_n, zxi1, zuvel_n, & |
|---|
| 188 | & zyj_n, zyj1, zvvel_n, ll_bounced ) |
|---|
| 189 | |
|---|
| 190 | pt%uvel = zuvel_n !** save in berg structure |
|---|
| 191 | pt%vvel = zvvel_n |
|---|
| 192 | pt%xi = zxi_n |
|---|
| 193 | pt%yj = zyj_n |
|---|
| 194 | |
|---|
| 195 | ! update actual position |
|---|
| 196 | pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) |
|---|
| 197 | pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) |
|---|
| 198 | |
|---|
| 199 | berg => berg%next ! switch to the next berg |
|---|
| 200 | ! |
|---|
| 201 | END DO !== end loop over all bergs ==! |
|---|
| 202 | ! |
|---|
| 203 | END SUBROUTINE icb_dyn |
|---|
| 204 | |
|---|
| 205 | |
|---|
| 206 | SUBROUTINE icb_ground( pi, pi0, pu, & |
|---|
| 207 | & pj, pj0, pv, ld_bounced ) |
|---|
| 208 | !!---------------------------------------------------------------------- |
|---|
| 209 | !! *** ROUTINE icb_ground *** |
|---|
| 210 | !! |
|---|
| 211 | !! ** Purpose : iceberg grounding. |
|---|
| 212 | !! |
|---|
| 213 | !! ** Method : - adjust velocity and then put iceberg back to start position |
|---|
| 214 | !! NB two possibilities available one of which is hard-coded here |
|---|
| 215 | !!---------------------------------------------------------------------- |
|---|
| 216 | REAL(wp), INTENT(inout) :: pi , pj ! current iceberg position |
|---|
| 217 | REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position |
|---|
| 218 | REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities |
|---|
| 219 | LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator |
|---|
| 220 | ! |
|---|
| 221 | INTEGER :: ii, ii0 |
|---|
| 222 | INTEGER :: ij, ij0 |
|---|
| 223 | INTEGER :: ibounce_method |
|---|
| 224 | !!---------------------------------------------------------------------- |
|---|
| 225 | ! |
|---|
| 226 | ld_bounced = .FALSE. |
|---|
| 227 | ! |
|---|
| 228 | ii0 = INT( pi0+0.5 ) ; ij0 = INT( pj0+0.5 ) ! initial gridpoint position (T-cell) |
|---|
| 229 | ii = INT( pi +0.5 ) ; ij = INT( pj +0.5 ) ! current - - |
|---|
| 230 | ! |
|---|
| 231 | IF( ii == ii0 .AND. ij == ij0 ) RETURN ! berg remains in the same cell |
|---|
| 232 | ! |
|---|
| 233 | ! map into current processor |
|---|
| 234 | ii0 = mi1( ii0 ) |
|---|
| 235 | ij0 = mj1( ij0 ) |
|---|
| 236 | ii = mi1( ii ) |
|---|
| 237 | ij = mj1( ij ) |
|---|
| 238 | ! |
|---|
| 239 | IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one |
|---|
| 240 | ! |
|---|
| 241 | ! From here, berg have reach land: treat grounding/bouncing |
|---|
| 242 | ! ------------------------------- |
|---|
| 243 | ld_bounced = .TRUE. |
|---|
| 244 | |
|---|
| 245 | !! not obvious what should happen now |
|---|
| 246 | !! if berg tries to enter a land box, the only location we can return it to is the start |
|---|
| 247 | !! position (pi0,pj0), since it has to be in a wet box to do any melting; |
|---|
| 248 | !! first option is simply to set whole velocity to zero and move back to start point |
|---|
| 249 | !! second option (suggested by gm) is only to set the velocity component in the (i,j) direction |
|---|
| 250 | !! of travel to zero; at a coastal boundary this has the effect of sliding the berg along the coast |
|---|
| 251 | |
|---|
| 252 | ibounce_method = 2 |
|---|
| 253 | SELECT CASE ( ibounce_method ) |
|---|
| 254 | CASE ( 1 ) |
|---|
| 255 | pi = pi0 |
|---|
| 256 | pj = pj0 |
|---|
| 257 | pu = 0._wp |
|---|
| 258 | pv = 0._wp |
|---|
| 259 | CASE ( 2 ) |
|---|
| 260 | IF( ii0 /= ii ) THEN |
|---|
| 261 | pi = pi0 ! return back to the initial position |
|---|
| 262 | pu = 0._wp ! zeroing of velocity in the direction of the grounding |
|---|
| 263 | ENDIF |
|---|
| 264 | IF( ij0 /= ij ) THEN |
|---|
| 265 | pj = pj0 ! return back to the initial position |
|---|
| 266 | pv = 0._wp ! zeroing of velocity in the direction of the grounding |
|---|
| 267 | ENDIF |
|---|
| 268 | END SELECT |
|---|
| 269 | ! |
|---|
| 270 | END SUBROUTINE icb_ground |
|---|
| 271 | |
|---|
| 272 | |
|---|
| 273 | SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax, & |
|---|
| 274 | & pyj, pe2, pvvel, pvvel0, pay, pdt ) |
|---|
| 275 | !!---------------------------------------------------------------------- |
|---|
| 276 | !! *** ROUTINE icb_accel *** |
|---|
| 277 | !! |
|---|
| 278 | !! ** Purpose : compute the iceberg acceleration. |
|---|
| 279 | !! |
|---|
| 280 | !! ** Method : - sum the terms in the momentum budget |
|---|
| 281 | !!---------------------------------------------------------------------- |
|---|
| 282 | TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg |
|---|
| 283 | REAL(wp) , INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential |
|---|
| 284 | REAL(wp) , INTENT(in ) :: puvel , pvvel ! berg velocity [m/s] |
|---|
| 285 | REAL(wp) , INTENT(in ) :: puvel0, pvvel0 ! initial berg velocity [m/s] |
|---|
| 286 | REAL(wp) , INTENT( out) :: pe1, pe2 ! horizontal scale factor at (xi,yj) |
|---|
| 287 | REAL(wp) , INTENT(inout) :: pax, pay ! berg acceleration |
|---|
| 288 | REAL(wp) , INTENT(in ) :: pdt ! berg time step |
|---|
| 289 | ! |
|---|
| 290 | REAL(wp), PARAMETER :: pp_alpha = 0._wp ! |
|---|
| 291 | REAL(wp), PARAMETER :: pp_beta = 1._wp ! |
|---|
| 292 | REAL(wp), PARAMETER :: pp_vel_lim =15._wp ! max allowed berg speed |
|---|
| 293 | REAL(wp), PARAMETER :: pp_accel_lim = 1.e-2_wp ! max allowed berg acceleration |
|---|
| 294 | REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! |
|---|
| 295 | ! |
|---|
| 296 | INTEGER :: itloop |
|---|
| 297 | REAL(wp) :: zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss |
|---|
| 298 | REAL(wp) :: zvo, zvi, zva, zvwave, zssh_y |
|---|
| 299 | REAL(wp) :: zff, zT, zD, zW, zL, zM, zF |
|---|
| 300 | REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad |
|---|
| 301 | REAL(wp) :: z_ocn, z_atm, z_ice |
|---|
| 302 | REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop |
|---|
| 303 | REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi |
|---|
| 304 | REAL(wp) :: zuveln, zvveln, zus, zvs |
|---|
| 305 | !!---------------------------------------------------------------------- |
|---|
| 306 | |
|---|
| 307 | ! Interpolate gridded fields to berg |
|---|
| 308 | nknberg = berg%number(1) |
|---|
| 309 | CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x, & |
|---|
| 310 | & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss ) |
|---|
| 311 | |
|---|
| 312 | zM = berg%current_point%mass |
|---|
| 313 | zT = berg%current_point%thickness ! total thickness |
|---|
| 314 | zD = ( rn_rho_bergs / pp_rho_seawater ) * zT ! draught (keel depth) |
|---|
| 315 | zF = zT - zD ! freeboard |
|---|
| 316 | zW = berg%current_point%width |
|---|
| 317 | zL = berg%current_point%length |
|---|
| 318 | |
|---|
| 319 | zhi = MIN( zhi , zD ) |
|---|
| 320 | zD_hi = MAX( 0._wp, zD-zhi ) |
|---|
| 321 | |
|---|
| 322 | ! Wave radiation |
|---|
| 323 | zuwave = zua - zuo ; zvwave = zva - zvo ! Use wind speed rel. to ocean for wave model |
|---|
| 324 | zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; |
|---|
| 325 | ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. |
|---|
| 326 | zampl = 0.5 * 0.02025 * zwmod ! This is "a", the wave amplitude |
|---|
| 327 | zLwavelength = 0.32 * zwmod ! Surface wave length fitted to data in table at |
|---|
| 328 | ! ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html |
|---|
| 329 | zLcutoff = 0.125 * zLwavelength |
|---|
| 330 | zLtop = 0.25 * zLwavelength |
|---|
| 331 | zCr = pp_Cr0 * MIN( MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.) ! Wave radiation coefficient |
|---|
| 332 | ! ! fitted to graph from Carrieres et al., POAC Drift Model. |
|---|
| 333 | zwave_rad = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) |
|---|
| 334 | zwmod = SQRT( zua*zua + zva*zva ) ! Wind speed |
|---|
| 335 | IF( zwmod /= 0._wp ) THEN |
|---|
| 336 | zuwave = zua/zwmod ! Wave radiation force acts in wind direction ... !!gm this should be the wind rel. to ocean ? |
|---|
| 337 | zvwave = zva/zwmod |
|---|
| 338 | ELSE |
|---|
| 339 | zuwave = 0. ; zvwave=0. ; zwave_rad=0. ! ... and only when wind is present. !!gm wave_rad=0. is useless |
|---|
| 340 | ENDIF |
|---|
| 341 | |
|---|
| 342 | ! Weighted drag coefficients |
|---|
| 343 | z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) |
|---|
| 344 | z_atm = pp_rho_air / zM * (0.5*pp_Cd_av*zW*zF +pp_Cd_ah*zW*zL) |
|---|
| 345 | z_ice = pp_rho_ice / zM * (0.5*pp_Cd_iv*zW*zhi ) |
|---|
| 346 | IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp |
|---|
| 347 | |
|---|
| 348 | zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel |
|---|
| 349 | ! |
|---|
| 350 | DO itloop = 1, 2 ! Iterate on drag coefficients |
|---|
| 351 | ! |
|---|
| 352 | zus = 0.5 * ( zuveln + puvel ) |
|---|
| 353 | zvs = 0.5 * ( zvveln + pvvel ) |
|---|
| 354 | zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) |
|---|
| 355 | zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) |
|---|
| 356 | zdrag_ice = z_ice * SQRT( (zus-zui)*(zus-zui) + (zvs-zvi)*(zvs-zvi) ) |
|---|
| 357 | ! |
|---|
| 358 | ! Explicit accelerations |
|---|
| 359 | !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & |
|---|
| 360 | ! -zdrag_ocn*(puvel-zuo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) |
|---|
| 361 | !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & |
|---|
| 362 | ! -zdrag_ocn*(pvvel-zvo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) |
|---|
| 363 | zaxe = -grav * zssh_x + zwave_rad * zuwave |
|---|
| 364 | zaye = -grav * zssh_y + zwave_rad * zvwave |
|---|
| 365 | IF( pp_alpha > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
|---|
| 366 | zaxe = zaxe + zff*pvvel0 |
|---|
| 367 | zaye = zaye - zff*puvel0 |
|---|
| 368 | ELSE |
|---|
| 369 | zaxe = zaxe + zff*pvvel |
|---|
| 370 | zaye = zaye - zff*puvel |
|---|
| 371 | ENDIF |
|---|
| 372 | IF( pp_beta > 0._wp ) THEN ! If implicit, use time-level (n) rather than RK4 latest |
|---|
| 373 | zaxe = zaxe - zdrag_ocn*(puvel0-zuo) - zdrag_atm*(puvel0-zua) -zdrag_ice*(puvel0-zui) |
|---|
| 374 | zaye = zaye - zdrag_ocn*(pvvel0-zvo) - zdrag_atm*(pvvel0-zva) -zdrag_ice*(pvvel0-zvi) |
|---|
| 375 | ELSE |
|---|
| 376 | zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui) |
|---|
| 377 | zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi) |
|---|
| 378 | ENDIF |
|---|
| 379 | |
|---|
| 380 | ! Solve for implicit accelerations |
|---|
| 381 | IF( pp_alpha + pp_beta > 0._wp ) THEN |
|---|
| 382 | zlambda = zdrag_ocn + zdrag_atm + zdrag_ice |
|---|
| 383 | zA11 = 1._wp + pp_beta *pdt*zlambda |
|---|
| 384 | zA12 = pp_alpha*pdt*zff |
|---|
| 385 | zdetA = 1._wp / ( zA11*zA11 + zA12*zA12 ) |
|---|
| 386 | pax = zdetA * ( zA11*zaxe + zA12*zaye ) |
|---|
| 387 | pay = zdetA * ( zA11*zaye - zA12*zaxe ) |
|---|
| 388 | ELSE |
|---|
| 389 | pax = zaxe ; pay = zaye |
|---|
| 390 | ENDIF |
|---|
| 391 | |
|---|
| 392 | zuveln = puvel0 + pdt*pax |
|---|
| 393 | zvveln = pvvel0 + pdt*pay |
|---|
| 394 | ! |
|---|
| 395 | END DO ! itloop |
|---|
| 396 | ! ! check the speed and acceleration limits |
|---|
| 397 | IF (nn_verbose_level > 0) THEN |
|---|
| 398 | IF( ABS( zuveln ) > pp_vel_lim .OR. ABS( zvveln ) > pp_vel_lim ) & |
|---|
| 399 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive velocity' |
|---|
| 400 | IF( ABS( pax ) > pp_accel_lim .OR. ABS( pay ) > pp_accel_lim ) & |
|---|
| 401 | WRITE(numicb,'("pe=",i3,x,a)') narea,'Dump triggered by excessive acceleration' |
|---|
| 402 | ENDIF |
|---|
| 403 | ! |
|---|
| 404 | END SUBROUTINE icb_accel |
|---|
| 405 | |
|---|
| 406 | SUBROUTINE icb_speed_limit( pe1, pe2, pdt, pu, pv, pn_speed_limit, pl_loop, pvel_factor, pn_count, pn_stage ) |
|---|
| 407 | !!---------------------------------------------------------------------- |
|---|
| 408 | !! *** ROUTINE icb_speed_limit *** |
|---|
| 409 | !! |
|---|
| 410 | !! ** Purpose : calculate CFL-based speed limit for icebergs. |
|---|
| 411 | !! |
|---|
| 412 | !! ** Method : - if iceberg velocity exceeds CFL-based limit, calculate |
|---|
| 413 | !! a reduction factor for the velocity. |
|---|
| 414 | !!---------------------------------------------------------------------- |
|---|
| 415 | REAL(wp), INTENT(in ) :: pe1, pe2 ! scale factors for this iceberg grid cell |
|---|
| 416 | REAL(wp), INTENT(in ) :: pdt ! timestep (for this RK4 stage) |
|---|
| 417 | REAL(wp), INTENT(inout) :: pu , pv ! current iceberg velocities |
|---|
| 418 | REAL(wp), INTENT(in ) :: pn_speed_limit ! CFL-based speed limit |
|---|
| 419 | LOGICAL , INTENT(inout) :: pl_loop ! flag to trigger recalculation of RK4 step |
|---|
| 420 | REAL(wp), INTENT(inout) :: pvel_factor ! factor to multiply velocities |
|---|
| 421 | INTEGER, INTENT(in ) :: pn_count ! speed limiting iteration that we are on |
|---|
| 422 | INTEGER, INTENT(in ) :: pn_stage ! RK4 stage that we are on |
|---|
| 423 | ! |
|---|
| 424 | REAL :: zspeed, zspeed_new, zloc_dx |
|---|
| 425 | !!---------------------------------------------------------------------- |
|---|
| 426 | |
|---|
| 427 | zspeed = SQRT( pu*pu + pv*pv ) ! Speed of berg |
|---|
| 428 | IF( zspeed > 0._wp ) THEN |
|---|
| 429 | zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing |
|---|
| 430 | zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt |
|---|
| 431 | IF( zspeed_new < zspeed ) THEN |
|---|
| 432 | pl_loop = .true. |
|---|
| 433 | ! arbitrary 0.99 factor to speed up convergence |
|---|
| 434 | pvel_factor = pvel_factor * 0.99_wp * zspeed_new/zspeed |
|---|
| 435 | CALL icb_dia_speed(pvel_factor, pn_count, pn_stage) |
|---|
| 436 | ENDIF |
|---|
| 437 | ENDIF |
|---|
| 438 | END SUBROUTINE icb_speed_limit |
|---|
| 439 | |
|---|
| 440 | !!====================================================================== |
|---|
| 441 | END MODULE icbdyn |
|---|
