Changeset 10364 for NEMO/trunk/src/OCE/DYN/sshwzv.F90
- Timestamp:
- 2018-11-30T18:42:51+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/DYN/sshwzv.F90
r10068 r10364 28 28 #endif 29 29 ! 30 USE iom 30 31 USE in_out_manager ! I/O manager 31 32 USE restart ! only for lrst_oce … … 41 42 PUBLIC ssh_nxt ! called by step.F90 42 43 PUBLIC wzv ! called by step.F90 44 PUBLIC wAimp ! called by step.F90 43 45 PUBLIC ssh_swp ! called by step.F90 44 46 … … 265 267 END SUBROUTINE ssh_swp 266 268 269 SUBROUTINE wAimp( kt ) 270 !!---------------------------------------------------------------------- 271 !! *** ROUTINE wAimp *** 272 !! 273 !! ** Purpose : compute the Courant number and partition vertical velocity 274 !! if a proportion needs to be treated implicitly 275 !! 276 !! ** Method : - 277 !! 278 !! ** action : wn : now vertical velocity (to be handled explicitly) 279 !! : wi : now vertical velocity (for implicit treatment) 280 !! 281 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 282 !!---------------------------------------------------------------------- 283 INTEGER, INTENT(in) :: kt ! time step 284 ! 285 INTEGER :: ji, jj, jk ! dummy loop indices 286 REAL(wp) :: zCu, zcff, z1_e3w ! local scalars 287 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 288 REAL(wp) , PARAMETER :: Cu_max = 0.27 ! local parameters 289 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 290 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 291 !!---------------------------------------------------------------------- 292 ! 293 IF( ln_timing ) CALL timing_start('wAimp') 294 ! 295 IF( kt == nit000 ) THEN 296 IF(lwp) WRITE(numout,*) 297 IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 298 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_n(ji,jj,jk) 307 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 308 & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 309 & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 310 & * r1_e1e2t(ji,jj) & 311 & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 312 & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 313 & * r1_e1e2t(ji,jj) & 314 & ) * z1_e3w 315 END DO 316 END DO 317 END DO 318 ! 319 CALL iom_put("Courant",Cu_adv) 320 ! 321 wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero 322 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 323 DO jk = 1, jpkm1 ! or scan Courant criterion and partition 324 DO jj = 2, jpjm1 ! w where necessary 325 DO ji = 2, fs_jpim1 ! vector opt. 326 ! 327 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 328 ! 329 IF( zCu < Cu_min ) THEN !<-- Fully explicit 330 zcff = 0._wp 331 ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit 332 zcff = ( zCu - Cu_min )**2 333 zcff = zcff / ( Fcu + zcff ) 334 ELSE !<-- Mostly implicit 335 zcff = ( zCu - Cu_max )/ zCu 336 ENDIF 337 zcff = MIN(1._wp, zcff) 338 ! 339 wi(ji,jj,jk) = zcff * wn(ji,jj,jk) 340 wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 341 ! 342 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient 343 END DO 344 END DO 345 END DO 346 ELSE 347 ! Fully explicit everywhere 348 Cu_adv = 0.0_wp ! Reuse array to output coefficient 349 ENDIF 350 CALL iom_put("wimp",wi) 351 CALL iom_put("wi_cff",Cu_adv) 352 CALL iom_put("wexp",wn) 353 ! 354 IF( ln_timing ) CALL timing_stop('wAimp') 355 ! 356 END SUBROUTINE wAimp 267 357 !!====================================================================== 268 358 END MODULE sshwzv
Note: See TracChangeset
for help on using the changeset viewer.