Ignore:
Timestamp:
2018-07-19T17:58:12+02:00 (3 years ago)
Author:
acc
Message:

Branch: dev_r9956_ENHANCE05_ZAD_AIMP. First set of changes to implement an adaptive-implicit vertical advection option (see ticket #2042). This code compiles and runs but has issues when the new option is activated.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90

    r9957 r9976  
    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)             ::   zcff_min, zcff_max                    ! local scalars 
     288      REAL(wp) , PARAMETER ::   Cu_min = 0.6_wp                       ! local parameters 
     289      REAL(wp) , PARAMETER ::   Cu_max = 0.95_wp                      ! local parameters 
     290      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
     291      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     292      !!---------------------------------------------------------------------- 
     293      ! 
     294      IF( ln_timing )   CALL timing_start('wAimp') 
     295      ! 
     296      IF( kt == nit000 ) THEN 
     297         IF(lwp) WRITE(numout,*) 
     298         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 
     299         IF(lwp) WRITE(numout,*) '~~~~~ ' 
     300         ! 
     301         wi    (:,:,jpk) = 0._wp              ! bottom value : wi=0 (set once for all) 
     302         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all) 
     303      ENDIF 
     304      ! 
     305      DO jk = 1, jpkm1            ! calculate Courant numbers 
     306         DO jj = 2, jpjm1 
     307            DO ji = 2, fs_jpim1   ! vector opt. 
     308               z1_e3w = 1._wp / e3w_n(ji,jj,jk) 
     309               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    & 
     310                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     311                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     312                  &                        / e2t(ji,jj)                                                       & 
     313                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     314                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     315                  &                        / e1t(ji,jj)                                                       & 
     316                  &                      ) * z1_e3w 
     317            END DO 
     318         END DO 
     319      END DO 
     320      ! 
     321      CALL iom_put("Courant",Cu_adv) 
     322      ! 
     323      wi(:,:,:) = 0._wp 
     324      zcff_min=999. ; zcff_max = -999. 
     325      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
     326         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition 
     327            DO jj = 2, jpjm1                            ! w where necessary 
     328               DO ji = 2, fs_jpim1   ! vector opt. 
     329                  ! 
     330                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 
     331                  ! 
     332                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit 
     333                     zcff = 0._wp 
     334                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
     335                     zcff = ( zCu - Cu_min )**2 
     336                     zcff = zcff / ( Fcu + zcff ) 
     337                  ELSE                                  !<-- Fully implicit 
     338                     zcff = ( zCu - Cu_max )/ zCu 
     339                  ENDIF 
     340                  ! 
     341                  !zcff = 0._wp 
     342                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk) 
     343                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
     344                  zcff_min = MIN( zcff, zcff_min ) 
     345                  zcff_max = MAX( zcff, zcff_max ) 
     346                  ! 
     347               END DO 
     348            END DO 
     349         END DO 
     350      ELSE 
     351         ! Fully explicit everywhere 
     352         write(*,*) 'NoWi ',kt 
     353      ENDIF 
     354         write(*,*) 'ZCFF ',kt,zcff_min,zcff_max 
     355      CALL iom_put("wimp",wi) 
     356      CALL iom_put("wexp",wn) 
     357      ! 
     358      IF( ln_timing )   CALL timing_stop('wAimp') 
     359      ! 
     360   END SUBROUTINE wAimp 
    267361   !!====================================================================== 
    268362END MODULE sshwzv 
Note: See TracChangeset for help on using the changeset viewer.