New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12149 for NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DYN/sshwzv.F90 – NEMO

Ignore:
Timestamp:
2019-12-10T15:03:24+01:00 (4 years ago)
Author:
ayoung
Message:

Updated trunk to 12072

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/ENHANCE-03_closea/src/OCE/DYN/sshwzv.F90

    r10907 r12149  
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1010   !!            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 
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    279280      !!            :   wi      : now vertical velocity (for implicit treatment) 
    280281      !! 
    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. 
    282285      !!---------------------------------------------------------------------- 
    283286      INTEGER, INTENT(in) ::   kt   ! time step 
    284287      ! 
    285288      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    286       REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars 
     289      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars 
    287290      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
    288       REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     291      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters 
    289292      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
    290293      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     
    297300         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
    298301         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 
    314321            END DO 
    315322         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 
    317341      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 
    318342      ! 
     
    320344      ! 
    321345      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
    322          DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
     346         DO jk = jpkm1, 2, -1                           ! or scan Courant criterion and partition 
    323347            DO jj = 1, jpj                              ! w where necessary 
    324348               DO ji = 1, jpi 
    325349                  ! 
    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  
    327357                  ! 
    328358                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
     
    339369                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
    340370                  ! 
    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 
    342372               END DO 
    343373            END DO 
    344374         END DO 
     375         Cu_adv(:,:,1) = 0._wp  
    345376      ELSE 
    346377         ! 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 
    348379         wi    (:,:,:) = 0._wp 
    349380      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.