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