Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
r11480 r11822 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 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 12 13 !!---------------------------------------------------------------------- … … 278 279 !! : wi : now vertical velocity (for implicit treatment) 279 280 !! 280 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 281 !! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 282 !! implicit scheme for vertical advection in oceanic modeling. 283 !! Ocean Modelling, 91, 38-69. 281 284 !!---------------------------------------------------------------------- 282 285 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 ! 300 Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all) 301 ENDIF 302 ! 303 DO jk = 1, jpkm1 ! calculate Courant numbers 304 DO jj = 2, jpjm1 305 DO ji = 2, fs_jpim1 ! vector opt. 306 z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm) 307 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 308 & + ( MAX( e2u(ji ,jj)*e3uw(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 309 & MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 310 & * r1_e1e2t(ji,jj) & 311 & + ( MAX( e1v(ji,jj )*e3vw(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 312 & MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 313 & * r1_e1e2t(ji,jj) & 314 & ) * 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(ji,jj,jk,Kmm) 311 ! 2*rdt and not r2dt (for restartability) 312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 313 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 314 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 315 & * r1_e1e2t(ji,jj) & 316 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 317 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 318 & * r1_e1e2t(ji,jj) & 319 & ) * z1_e3t 320 END DO 315 321 END DO 316 322 END DO 317 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(ji,jj,jk,Kmm) 328 ! 2*rdt and not r2dt (for restartability) 329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 330 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 331 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 332 & * r1_e1e2t(ji,jj) & 333 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 334 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 335 & * r1_e1e2t(ji,jj) & 336 & ) * z1_e3t 337 END DO 338 END DO 339 END DO 340 ENDIF 341 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 318 342 ! 319 343 CALL iom_put("Courant",Cu_adv) 320 344 ! 321 wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero322 345 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 323 DO jk = 1, jpkm1! or scan Courant criterion and partition324 DO jj = 2, jpjm1! w where necessary325 DO ji = 2, fs_jpim1 ! vector opt.346 DO jk = jpkm1, 2, -1 ! or scan Courant criterion and partition 347 DO jj = 1, jpj ! w where necessary 348 DO ji = 1, jpi 326 349 ! 327 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 328 357 ! 329 IF( zCu < Cu_min ) THEN!<-- Fully explicit358 IF( zCu <= Cu_min ) THEN !<-- Fully explicit 330 359 zcff = 0._wp 331 360 ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit … … 340 369 ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 341 370 ! 342 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 343 372 END DO 344 373 END DO 345 374 END DO 375 Cu_adv(:,:,1) = 0._wp 346 376 ELSE 347 377 ! Fully explicit everywhere 348 Cu_adv = 0.0_wp ! Reuse array to output coefficient 378 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl 379 wi (:,:,:) = 0._wp 349 380 ENDIF 350 381 CALL iom_put("wimp",wi)
Note: See TracChangeset
for help on using the changeset viewer.
