- Timestamp:
- 2020-12-14T15:50:33+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdyn.F90
r14102 r14167 41 41 INTEGER, INTENT(in) :: kt ! 42 42 ! 43 LOGICAL :: ll_bounced 43 LOGICAL :: ll_bounced, l_loop 44 44 REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 45 45 REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 … … 47 47 REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 48 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) 49 51 REAL(wp) :: zspeed, zspeed_new, zloc_dx 50 52 REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 53 INTEGER :: icount 51 54 TYPE(iceberg), POINTER :: berg 52 55 TYPE(point) , POINTER :: pt … … 77 80 ll_bounced = .FALSE. 78 81 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( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & 88 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 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( zxi2, zxi1, zu1, & 101 & zyj2, zyj1, zv1, ll_bounced ) 102 103 ! !** A2 = A(X2,V2) 104 CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & 105 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 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( zxi3, zxi1, zu3, & 117 & zyj3, zyj1, zv3, ll_bounced ) 118 119 ! !** A3 = A(X3,V3) 120 CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & 121 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 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( zxi4, zxi1, zu4, & 133 & zyj4, zyj1, zv4, ll_bounced ) 134 135 ! !** A4 = A(X4,V4) 136 CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & 137 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 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 ) 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) 150 186 151 187 CALL icb_ground( zxi_n, zxi1, zuvel_n, & 152 188 & zyj_n, zyj1, zvvel_n, ll_bounced ) 153 154 IF( rn_speed_limit > 0._wp ) THEN ! Limit speed of bergs based on a CFL criteria (if asked)155 zspeed = SQRT( zuvel_n*zuvel_n + zvvel_n*zvvel_n ) ! Speed of berg156 IF( zspeed > 0._wp ) THEN157 zloc_dx = MIN( ze1, ze2 ) ! minimum grid spacing158 zspeed_new = zloc_dx / zdt * rn_speed_limit ! Speed limit as a factor of dx / dt159 IF( zspeed_new < zspeed ) THEN160 zuvel_n = zuvel_n * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed161 zvvel_n = zvvel_n * ( zspeed_new / zspeed ) ! without changing the direction162 CALL icb_dia_speed()163 ENDIF164 ENDIF165 ENDIF166 189 167 190 pt%uvel = zuvel_n !** save in berg structure … … 381 404 END SUBROUTINE icb_accel 382 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 383 440 !!====================================================================== 384 441 END MODULE icbdyn
Note: See TracChangeset
for help on using the changeset viewer.