- Timestamp:
- 2019-12-10T15:03:24+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DYN/sshwzv.F90
r10907 r12149 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection 11 12 !!---------------------------------------------------------------------- 12 13 … … 279 280 !! : wi : now vertical velocity (for implicit treatment) 280 281 !! 281 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 282 !! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 283 !! implicit scheme for vertical advection in oceanic modeling. 284 !! Ocean Modelling, 91, 38-69. 282 285 !!---------------------------------------------------------------------- 283 286 INTEGER, INTENT(in) :: kt ! time step 284 287 ! 285 288 INTEGER :: ji, jj, jk ! dummy loop indices 286 REAL(wp) :: zCu, zcff, z1_e3 w! local scalars289 REAL(wp) :: zCu, zcff, z1_e3t ! local scalars 287 290 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 288 REAL(wp) , PARAMETER :: Cu_max = 0. 27! local parameters291 REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters 289 292 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 290 293 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters … … 297 300 IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 298 301 IF(lwp) WRITE(numout,*) '~~~~~ ' 299 ENDIF 300 ! 301 ! 302 DO jk = 1, jpkm1 ! calculate Courant numbers 303 DO jj = 2, jpjm1 304 DO ji = 2, fs_jpim1 ! vector opt. 305 z1_e3w = 1._wp / e3w_n(ji,jj,jk) 306 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & ! 2*rdt and not r2dt (for restartability) 307 & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 308 & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 309 & * r1_e1e2t(ji,jj) & 310 & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 311 & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 312 & * r1_e1e2t(ji,jj) & 313 & ) * z1_e3w 302 wi(:,:,:) = 0._wp 303 ENDIF 304 ! 305 ! Calculate Courant numbers 306 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 307 DO jk = 1, jpkm1 308 DO jj = 2, jpjm1 309 DO ji = 2, fs_jpim1 ! vector opt. 310 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 311 ! 2*rdt and not r2dt (for restartability) 312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 313 & + ( MAX( e2u(ji ,jj)*e3u_n(ji ,jj,jk)*un(ji ,jj,jk) + un_td(ji ,jj,jk), 0._wp ) - & 314 & MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) ) & 315 & * r1_e1e2t(ji,jj) & 316 & + ( MAX( e1v(ji,jj )*e3v_n(ji,jj ,jk)*vn(ji,jj ,jk) + vn_td(ji,jj ,jk), 0._wp ) - & 317 & MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) ) & 318 & * r1_e1e2t(ji,jj) & 319 & ) * z1_e3t 320 END DO 314 321 END DO 315 322 END DO 316 END DO 323 ELSE 324 DO jk = 1, jpkm1 325 DO jj = 2, jpjm1 326 DO ji = 2, fs_jpim1 ! vector opt. 327 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 328 ! 2*rdt and not r2dt (for restartability) 329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 330 & + ( MAX( e2u(ji ,jj)*e3u_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 331 & MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 332 & * r1_e1e2t(ji,jj) & 333 & + ( MAX( e1v(ji,jj )*e3v_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 334 & MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 335 & * r1_e1e2t(ji,jj) & 336 & ) * z1_e3t 337 END DO 338 END DO 339 END DO 340 ENDIF 317 341 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 318 342 ! … … 320 344 ! 321 345 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 322 DO jk = 1, jpkm1! or scan Courant criterion and partition346 DO jk = jpkm1, 2, -1 ! or scan Courant criterion and partition 323 347 DO jj = 1, jpj ! w where necessary 324 348 DO ji = 1, jpi 325 349 ! 326 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 350 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 351 ! alt: 352 ! IF ( wn(ji,jj,jk) > 0._wp ) THEN 353 ! zCu = Cu_adv(ji,jj,jk) 354 ! ELSE 355 ! zCu = Cu_adv(ji,jj,jk-1) 356 ! ENDIF 327 357 ! 328 358 IF( zCu <= Cu_min ) THEN !<-- Fully explicit … … 339 369 wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 340 370 ! 341 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient 371 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl 342 372 END DO 343 373 END DO 344 374 END DO 375 Cu_adv(:,:,1) = 0._wp 345 376 ELSE 346 377 ! Fully explicit everywhere 347 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient 378 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl 348 379 wi (:,:,:) = 0._wp 349 380 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.