Ignore:
Timestamp:
2018-11-30T18:42:51+01:00 (2 years ago)
Author:
acc
Message:

Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.

File:
1 edited

Legend:

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

    r10068 r10364  
    2828#endif 
    2929   ! 
     30   USE iom  
    3031   USE in_out_manager ! I/O manager 
    3132   USE restart        ! only for lrst_oce 
     
    4142   PUBLIC   ssh_nxt    ! called by step.F90 
    4243   PUBLIC   wzv        ! called by step.F90 
     44   PUBLIC   wAimp      ! called by step.F90 
    4345   PUBLIC   ssh_swp    ! called by step.F90 
    4446 
     
    265267   END SUBROUTINE ssh_swp 
    266268 
     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 
    267357   !!====================================================================== 
    268358END MODULE sshwzv 
Note: See TracChangeset for help on using the changeset viewer.