Changeset 11293


Ignore:
Timestamp:
2019-07-18T15:21:40+02:00 (13 months ago)
Author:
jchanut
Message:

#2292, Fixes courant number computation in adaptive vertical advection

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/DYN/sshwzv.F90

    r10907 r11293  
    284284      ! 
    285285      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    286       REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars 
     286      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars 
    287287      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
    288288      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     
    297297         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
    298298         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     299         wi(:,:,:) = 0._wp 
    299300      ENDIF 
    300301      ! 
     
    303304         DO jj = 2, jpjm1 
    304305            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 
     306               z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     307               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) 
     308                  &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     309                  &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     310                  &                               * r1_e1e2t(ji,jj)                                                 & 
     311                  &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     312                  &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     313                  &                               * r1_e1e2t(ji,jj)                                                 & 
     314                  &                             ) * z1_e3t 
    314315            END DO 
    315316         END DO 
     
    320321      ! 
    321322      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
    322          DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
     323         DO jk = jpkm1, 2, -1                           ! or scan Courant criterion and partition 
    323324            DO jj = 1, jpj                              ! w where necessary 
    324325               DO ji = 1, jpi 
    325326                  ! 
    326                   zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 
     327                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 
     328! alt: 
     329!                  IF ( wn(ji,jj,jk) > 0._wp ) THEN  
     330!                     zCu =  Cu_adv(ji,jj,jk)  
     331!                  ELSE 
     332!                     zCu =  Cu_adv(ji,jj,jk-1) 
     333!                  ENDIF  
    327334                  ! 
    328335                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
     
    343350            END DO 
    344351         END DO 
     352         Cu_adv(:,:,1) = 0._wp  
    345353      ELSE 
    346354         ! Fully explicit everywhere 
Note: See TracChangeset for help on using the changeset viewer.