Changeset 14274
- Timestamp:
- 2021-01-06T12:23:06+01:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdia.F90
r13587 r14274 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 89 INTEGER , DIMENSION(4) :: nspeeding_tickets 89 90 INTEGER , DIMENSION(nclasses) :: nbergs_calved_by_class 90 91 … … 124 125 nbergs_calved = 0 125 126 nbergs_calved_by_class(:) = 0 126 nspeeding_tickets = 0 127 nspeeding_tickets(:) = 0 128 vel_factor_min = 1._wp 127 129 stored_heat_end = 0._wp 128 130 floating_heat_end = 0._wp … … 157 159 IF( lk_mpp ) THEN 158 160 ALLOCATE( rsumbuf(23) ) ; rsumbuf(:) = 0._wp 159 ALLOCATE( nsumbuf( 4+nclasses) ) ; nsumbuf(:) = 0161 ALLOCATE( nsumbuf(7+nclasses) ) ; nsumbuf(:) = 0 160 162 rsumbuf(1) = floating_mass_start 161 163 rsumbuf(2) = bergs_mass_start … … 265 267 nsumbuf(2) = nbergs_calved 266 268 nsumbuf(3) = nbergs_melted 267 nsumbuf(4 ) = nspeeding_tickets269 nsumbuf(4:7) = nspeeding_tickets(:) 268 270 DO ik = 1, nclasses 269 nsumbuf( 4+ik) = nbergs_calved_by_class(ik)271 nsumbuf(7+ik) = nbergs_calved_by_class(ik) 270 272 END DO 271 CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+ 4), nclasses+4)273 CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+7), nclasses+7 ) 272 274 ! 273 275 nbergs_end = nsumbuf(1) 274 276 nbergs_calved = nsumbuf(2) 275 277 nbergs_melted = nsumbuf(3) 276 nspeeding_tickets = nsumbuf(4)278 nspeeding_tickets(:) = nsumbuf(4:7) 277 279 DO ik = 1,nclasses 278 nbergs_calved_by_class(ik)= nsumbuf( 4+ik)280 nbergs_calved_by_class(ik)= nsumbuf(7+ik) 279 281 END DO 280 282 ! 283 CALL mpp_min( 'icbdia', vel_factor_min, 1 ) 281 284 ENDIF 282 285 ! … … 329 332 IF (nn_verbose_level > 0) THEN 330 333 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 334 WRITE( numicb, '("n speeding tickets by RK4 stage = ",i6,3(",",i6))') (nspeeding_tickets(ik),ik=1,4) 335 IF( SUM(nspeeding_tickets) > 0 ) THEN 336 WRITE( numicb, '("min velocity reduction factor = ",f12.8)') vel_factor_min 337 ENDIF 332 338 ENDIF 333 339 ! … … 337 343 nbergs_calved = 0 338 344 nbergs_calved_by_class(:) = 0 339 nspeeding_tickets = 0 345 nspeeding_tickets(:) = 0 346 vel_factor_min = 1._wp 340 347 stored_heat_start = stored_heat_end 341 348 floating_heat_start = floating_heat_end … … 474 481 475 482 476 SUBROUTINE icb_dia_speed() 477 !!---------------------------------------------------------------------- 478 !!---------------------------------------------------------------------- 479 ! 480 IF( .NOT.ln_bergdia ) RETURN 481 nspeeding_tickets = nspeeding_tickets + 1 483 SUBROUTINE icb_dia_speed(pvel_factor, pn_stage) 484 !!---------------------------------------------------------------------- 485 !!---------------------------------------------------------------------- 486 REAL(wp), INTENT(in) :: pvel_factor ! factor by which velocity reduced 487 INTEGER , INTENT(in) :: pn_stage ! which stage of the RK4 calculation are we on 488 ! 489 IF( (.NOT.ln_bergdia) .OR. pn_stage .lt. 1 .OR. pn_stage .gt. 4 ) RETURN 490 nspeeding_tickets(pn_stage) = nspeeding_tickets(pn_stage) + 1 491 vel_factor_min = MIN(vel_factor_min,pvel_factor) ! keep track of minimum reduction factor 482 492 ! 483 493 END SUBROUTINE icb_dia_speed -
NEMO/branches/UKMO/NEMO_4.0.3_icb_speed_limit2/src/OCE/ICB/icbdyn.F90
r13587 r14274 85 85 ! !** A1 = A(X1,V1) 86 86 CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & 87 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )87 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 1 ) 88 88 ! 89 89 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt … … 102 102 ! !** A2 = A(X2,V2) 103 103 CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & 104 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )104 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 2 ) 105 105 ! 106 106 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt … … 118 118 ! !** A3 = A(X3,V3) 119 119 CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & 120 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt )120 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 3 ) 121 121 ! 122 122 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt … … 134 134 ! !** A4 = A(X4,V4) 135 135 CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & 136 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt )136 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 4 ) 137 137 138 138 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt … … 235 235 236 236 SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax, & 237 & pyj, pe2, pvvel, pvvel0, pay, pdt )237 & pyj, pe2, pvvel, pvvel0, pay, pdt, pstage ) 238 238 !!---------------------------------------------------------------------- 239 239 !! *** ROUTINE icb_accel *** … … 250 250 REAL(wp) , INTENT(inout) :: pax, pay ! berg acceleration 251 251 REAL(wp) , INTENT(in ) :: pdt ! berg time step 252 INTEGER , INTENT(in ) :: pstage ! RK4 stage (for speed limiter) 252 253 ! 253 254 REAL(wp), PARAMETER :: pp_alpha = 0._wp ! … … 265 266 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 266 267 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 267 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 268 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale 268 269 !!---------------------------------------------------------------------- 269 270 … … 361 362 zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg 362 363 IF( zspeed > 0._wp ) THEN 363 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 364 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt 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 365 372 IF( zspeed_new < zspeed ) THEN 366 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed 367 zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction 368 CALL icb_dia_speed() 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) 369 381 ENDIF 370 382 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.