Changeset 14167
- Timestamp:
- 2020-12-14T15:50:33+01:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB
- Files:
-
- 2 edited
-
icbdia.F90 (modified) (7 diffs)
-
icbdyn.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit/src/OCE/ICB/icbdia.F90
r13587 r14167 83 83 REAL(wp) :: heat_to_bergs_net, heat_to_ocean_net, melt_net 84 84 REAL(wp) :: calving_to_bergs_net 85 REAL(wp) :: vel_factor_min 85 86 86 87 INTEGER :: nbergs_start, nbergs_end, nbergs_calved 87 88 INTEGER :: nbergs_melted 88 INTEGER :: nspeeding_tickets 89 INTEGER , DIMENSION(4) :: nspeeding_tickets 90 INTEGER :: icount_max 89 91 INTEGER , DIMENSION(nclasses) :: nbergs_calved_by_class 90 92 … … 124 126 nbergs_calved = 0 125 127 nbergs_calved_by_class(:) = 0 126 nspeeding_tickets = 0 128 nspeeding_tickets(:) = 0 129 icount_max = 0 130 vel_factor_min = 1._wp 127 131 stored_heat_end = 0._wp 128 132 floating_heat_end = 0._wp … … 157 161 IF( lk_mpp ) THEN 158 162 ALLOCATE( rsumbuf(23) ) ; rsumbuf(:) = 0._wp 159 ALLOCATE( nsumbuf( 4+nclasses) ) ; nsumbuf(:) = 0163 ALLOCATE( nsumbuf(7+nclasses) ) ; nsumbuf(:) = 0 160 164 rsumbuf(1) = floating_mass_start 161 165 rsumbuf(2) = bergs_mass_start … … 265 269 nsumbuf(2) = nbergs_calved 266 270 nsumbuf(3) = nbergs_melted 267 nsumbuf(4 ) = nspeeding_tickets271 nsumbuf(4:7) = nspeeding_tickets(:) 268 272 DO ik = 1, nclasses 269 nsumbuf( 4+ik) = nbergs_calved_by_class(ik)273 nsumbuf(7+ik) = nbergs_calved_by_class(ik) 270 274 END DO 271 CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+ 4), nclasses+4)275 CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+7), nclasses+7 ) 272 276 ! 273 277 nbergs_end = nsumbuf(1) 274 278 nbergs_calved = nsumbuf(2) 275 279 nbergs_melted = nsumbuf(3) 276 nspeeding_tickets = nsumbuf(4)280 nspeeding_tickets(:) = nsumbuf(4:7) 277 281 DO ik = 1,nclasses 278 nbergs_calved_by_class(ik)= nsumbuf( 4+ik)282 nbergs_calved_by_class(ik)= nsumbuf(7+ik) 279 283 END DO 280 284 ! 285 CALL mpp_min( 'icbdia', vel_factor_min, 1 ) 286 CALL mpp_max( 'icbdia', icount_max, 1 ) 281 287 ENDIF 282 288 ! … … 329 335 IF (nn_verbose_level > 0) THEN 330 336 WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses) 331 IF( nspeeding_tickets > 0 ) WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets 337 WRITE( numicb, '("n speeding tickets by RK4 stage = ",i6,3(",",i6))') (nspeeding_tickets(ik),ik=1,4) 338 IF( SUM(nspeeding_tickets) > 0 ) THEN 339 WRITE( numicb, '("max number of iterations = ",i6)') icount_max 340 WRITE( numicb, '("min velocity reduction factor = ",f12.8)') vel_factor_min 341 ENDIF 332 342 ENDIF 333 343 ! … … 337 347 nbergs_calved = 0 338 348 nbergs_calved_by_class(:) = 0 339 nspeeding_tickets = 0 349 nspeeding_tickets(:) = 0 350 icount_max = 0 351 vel_factor_min = 1._wp 340 352 stored_heat_start = stored_heat_end 341 353 floating_heat_start = floating_heat_end … … 474 486 475 487 476 SUBROUTINE icb_dia_speed() 477 !!---------------------------------------------------------------------- 478 !!---------------------------------------------------------------------- 479 ! 480 IF( .NOT.ln_bergdia ) RETURN 481 nspeeding_tickets = nspeeding_tickets + 1 488 SUBROUTINE icb_dia_speed(pvel_factor, pn_count, pn_stage) 489 !!---------------------------------------------------------------------- 490 !!---------------------------------------------------------------------- 491 REAL(wp), INTENT(in) :: pvel_factor ! factor by which velocity reduced 492 INTEGER , INTENT(in) :: pn_count ! which velocity-limiting stage are we on 493 INTEGER , INTENT(in) :: pn_stage ! which stage of the RK4 calculation are we on 494 ! 495 IF( (.NOT.ln_bergdia) .OR. pn_stage .lt. 1 .OR. pn_stage .gt. 4 ) RETURN 496 nspeeding_tickets(pn_stage) = nspeeding_tickets(pn_stage) + 1 497 vel_factor_min = MIN(vel_factor_min,pvel_factor) ! keep track of minimum reduction factor 498 icount_max = MAX(icount_max,pn_count) ! keep track of maximum number of iterations 482 499 ! 483 500 END SUBROUTINE icb_dia_speed -
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.
