Changeset 4736


Ignore:
Timestamp:
2014-08-05T12:24:00+02:00 (6 years ago)
Author:
acc
Message:

Branch 2014/dev_r4743_NOC2_ZTS, #1367. Code changes to introduce optional sub-timestepping for vertical advection. Changes in dynadv.F90, dynzad.F90, traadv.F90, traadv_tvd.F90 and namelist_ref only.

Location:
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4699 r4736  
    684684   ln_traadv_qck    =  .false.  !  QUICKEST scheme 
    685685   ln_traadv_msc_ups=  .false.  !  use upstream scheme within muscl 
     686   ln_traadv_tvd_zts=  .false.  !  TVD scheme with sub-timestepping of vertical tracer advection 
    686687/ 
    687688!----------------------------------------------------------------------- 
     
    758759   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme 
    759760   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
     761   ln_dynzad_zts = .false. !  Use (T) sub timestepping for vertical momentum advection 
    760762/ 
    761763!----------------------------------------------------------------------- 
  • branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r4624 r4736  
    3030   LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag 
    3131   LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag 
     32   LOGICAL, PUBLIC ::   ln_dynzad_zts   !: vertical advection with sub-timestepping (requires vector form) 
    3233    
    3334   INTEGER ::   nadv   ! choice of the formulation and scheme for the advection 
     
    6465                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy 
    6566                      CALL dyn_zad     ( kt )    ! vector form : vertical advection 
    66       CASE ( 1 )  
     67      CASE ( 1 )      
     68                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy 
     69                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping 
     70      CASE ( 2 )  
    6771                      CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme 
    68       CASE ( 2 )    
     72      CASE ( 3 )    
    6973                      CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme 
    7074      ! 
     
    9195      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    9296      !! 
    93       NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs 
     97      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 
    9498      !!---------------------------------------------------------------------- 
    9599 
     
    111115         WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 
    112116         WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs 
     117         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts 
    113118      ENDIF 
    114119 
     
    120125 
    121126      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 
     127      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   & 
     128          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 
    122129 
    123130      !                               ! Set nadv 
    124131      IF( ln_dynadv_vec  )   nadv =  0  
    125       IF( ln_dynadv_cen2 )   nadv =  1 
    126       IF( ln_dynadv_ubs  )   nadv =  2 
     132      IF( ln_dynzad_zts  )   nadv =  1 
     133      IF( ln_dynadv_cen2 )   nadv =  2 
     134      IF( ln_dynadv_ubs  )   nadv =  3 
    127135      IF( lk_esopa       )   nadv = -1 
    128136 
     
    130138         WRITE(numout,*) 
    131139         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used' 
    132          IF( nadv ==  1 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used' 
    133          IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : UBS       scheme is used' 
     140         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used' 
     141         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used' 
     142         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used' 
    134143         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation' 
    135144      ENDIF 
  • branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r3294 r4736  
    2727   PRIVATE 
    2828    
    29    PUBLIC   dyn_zad   ! routine called by step.F90 
     29   PUBLIC   dyn_zad       ! routine called by dynadv.F90 
     30   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90 
    3031 
    3132   !! * Substitutions 
     
    8384         DO jj = 2, jpj                   ! vertical fluxes  
    8485            DO ji = fs_2, jpi             ! vector opt. 
    85                zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 
     86               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 
    8687            END DO 
    8788         END DO 
     
    9596      DO jj = 2, jpjm1              ! Surface and bottom values set to zero 
    9697         DO ji = fs_2, fs_jpim1           ! vector opt. 
    97             zwuw(ji,jj, 1 ) = 0.e0 
    98             zwvw(ji,jj, 1 ) = 0.e0 
    99             zwuw(ji,jj,jpk) = 0.e0 
    100             zwvw(ji,jj,jpk) = 0.e0 
     98            zwuw(ji,jj, 1 ) = 0._wp 
     99            zwvw(ji,jj, 1 ) = 0._wp 
     100            zwuw(ji,jj,jpk) = 0._wp 
     101            zwvw(ji,jj,jpk) = 0._wp 
    101102         END DO   
    102103      END DO 
     
    132133   END SUBROUTINE dyn_zad 
    133134 
     135   SUBROUTINE dyn_zad_zts ( kt ) 
     136      !!---------------------------------------------------------------------- 
     137      !!                  ***  ROUTINE dynzad_zts  *** 
     138      !!  
     139      !! ** Purpose :   Compute the now vertical momentum advection trend and  
     140      !!      add it to the general trend of momentum equation. This version 
     141      !!      uses sub-timesteps for improved numerical stability with small 
     142      !!      vertical grid sizes. This is especially relevant when using  
     143      !!      embedded ice with thin surface boxes. 
     144      !! 
     145      !! ** Method  :   The now vertical advection of momentum is given by: 
     146      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 
     147      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 
     148      !!      Add this trend to the general trend (ua,va): 
     149      !!         (ua,va) = (ua,va) + w dz(u,v) 
     150      !! 
     151      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends 
     152      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 
     153      !!---------------------------------------------------------------------- 
     154      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx 
     155      ! 
     156      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices 
     157      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection 
     158      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps 
     159      REAL(wp) ::   zua, zva        ! temporary scalars 
     160      REAL(wp) ::   zr_rdt          ! temporary scalar 
     161      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection 
     162      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection 
     163      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww 
     164      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv 
     165      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs 
     166      !!---------------------------------------------------------------------- 
     167      ! 
     168      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts') 
     169      ! 
     170      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww )  
     171      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs )  
     172      ! 
     173      IF( kt == nit000 ) THEN 
     174         IF(lwp)WRITE(numout,*) 
     175         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps' 
     176      ENDIF 
     177 
     178      IF( l_trddyn )   THEN         ! Save ua and va trends 
     179         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
     180         ztrdu(:,:,:) = ua(:,:,:)  
     181         ztrdv(:,:,:) = va(:,:,:)  
     182      ENDIF 
     183       
     184      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     185          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping) 
     186      ELSE 
     187          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog) 
     188      ENDIF 
     189       
     190      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes 
     191         DO jj = 2, jpj                    
     192            DO ji = fs_2, jpi             ! vector opt. 
     193               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 
     194            END DO 
     195         END DO 
     196      END DO 
     197 
     198      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero 
     199         DO ji = fs_2, fs_jpim1           ! vector opt. 
     200            zwuw(ji,jj, 1 ) = 0._wp 
     201            zwvw(ji,jj, 1 ) = 0._wp 
     202            zwuw(ji,jj,jpk) = 0._wp 
     203            zwvw(ji,jj,jpk) = 0._wp 
     204         END DO   
     205      END DO 
     206 
     207! Start with before values and use sub timestepping to reach after values 
     208 
     209      zus(:,:,:,1) = ub(:,:,:) 
     210      zvs(:,:,:,1) = vb(:,:,:) 
     211 
     212      DO jl = 1, jnzts                   ! Start of sub timestepping loop 
     213 
     214         IF( jl == 1 ) THEN              ! Euler forward to kick things off 
     215           jtb = 1   ;   jtn = 1   ;   jta = 2 
     216           zts = z2dtzts 
     217         ELSEIF( jl == 2 ) THEN          ! First leapfrog step 
     218           jtb = 1   ;   jtn = 2   ;   jta = 3 
     219           zts = 2._wp * z2dtzts 
     220         ELSE                            ! Shuffle pointers for subsequent leapfrog steps 
     221           jtb = MOD(jtb,3) + 1 
     222           jtn = MOD(jtn,3) + 1 
     223           jta = MOD(jta,3) + 1 
     224         ENDIF 
     225 
     226         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical 
     227            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
     228               DO ji = fs_2, fs_jpim1        ! vector opt. 
     229                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 
     230                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 
     231               END DO   
     232            END DO    
     233         END DO 
     234         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points 
     235            DO jj = 2, jpjm1 
     236               DO ji = fs_2, fs_jpim1       ! vector opt. 
     237                  !                         ! vertical momentum advective trends 
     238                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     239                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
     240                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 
     241                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts 
     242               END DO   
     243            END DO   
     244         END DO 
     245 
     246      END DO      ! End of sub timestepping loop 
     247 
     248      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts ) 
     249      DO jk = 1, jpkm1              ! Recover trends over the outer timestep 
     250         DO jj = 2, jpjm1 
     251            DO ji = fs_2, fs_jpim1       ! vector opt. 
     252               !                         ! vertical momentum advective trends 
     253               !                         ! add the trends to the general momentum trends 
     254               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt 
     255               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt 
     256            END DO   
     257         END DO   
     258      END DO 
     259 
     260      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic 
     261         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     262         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     263         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt) 
     264         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
     265      ENDIF 
     266      !                             ! Control print 
     267      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   & 
     268         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     269      ! 
     270      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww )  
     271      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs )  
     272      ! 
     273      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts') 
     274      ! 
     275   END SUBROUTINE dyn_zad_zts 
     276 
    134277   !!====================================================================== 
    135278END MODULE dynzad 
  • branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r4624 r4736  
    4343   LOGICAL ::   ln_traadv_cen2     ! 2nd order centered scheme flag 
    4444   LOGICAL ::   ln_traadv_tvd      ! TVD scheme flag 
     45   LOGICAL ::   ln_traadv_tvd_zts  ! TVD scheme flag with vertical sub time-stepping 
    4546   LOGICAL ::   ln_traadv_muscl    ! MUSCL scheme flag 
    4647   LOGICAL ::   ln_traadv_muscl2   ! MUSCL2 scheme flag 
     
    120121      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS  
    121122      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST  
     123      CASE ( 7 )   ;   CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS 
    122124      ! 
    123125      CASE (-1 )                                      !==  esopa: test all possibility with control print  ==! 
     
    166168         &                 ln_traadv_muscl, ln_traadv_muscl2,  & 
    167169         &                 ln_traadv_ubs  , ln_traadv_qck,     & 
    168          &                 ln_traadv_msc_ups 
     170         &                 ln_traadv_msc_ups, ln_traadv_tvd_zts 
    169171      !!---------------------------------------------------------------------- 
    170172 
     
    190192         WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck     = ', ln_traadv_qck 
    191193         WRITE(numout,*) '      upstream scheme within muscl   ln_traadv_msc_ups = ', ln_traadv_msc_ups 
     194         WRITE(numout,*) '      TVD advection scheme with zts  ln_traadv_tvd_zts = ', ln_traadv_tvd_zts 
    192195      ENDIF 
    193196 
     
    199202      IF( ln_traadv_ubs    )   ioptio = ioptio + 1 
    200203      IF( ln_traadv_qck    )   ioptio = ioptio + 1 
     204      IF( ln_traadv_tvd_zts)   ioptio = ioptio + 1 
    201205      IF( lk_esopa         )   ioptio =          1 
    202206 
     
    210214      IF( ln_traadv_ubs    )   nadv =  5 
    211215      IF( ln_traadv_qck    )   nadv =  6 
     216      IF( ln_traadv_tvd_zts)   nadv =  7 
    212217      IF( lk_esopa         )   nadv = -1 
    213218 
     
    220225         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used' 
    221226         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used' 
     227         IF( nadv ==  7 )   WRITE(numout,*) '         TVD ZTS   scheme is used' 
    222228         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme' 
    223229      ENDIF 
  • branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4499 r4736  
    3737   PRIVATE 
    3838 
    39    PUBLIC   tra_adv_tvd    ! routine called by step.F90 
     39   PUBLIC   tra_adv_tvd        ! routine called by traadv.F90 
     40   PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90 
    4041 
    4142   LOGICAL ::   l_trd   ! flag to compute trends 
     
    247248   END SUBROUTINE tra_adv_tvd 
    248249 
     250   SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
     251      &                                       ptb, ptn, pta, kjpt ) 
     252      !!---------------------------------------------------------------------- 
     253      !!                  ***  ROUTINE tra_adv_tvd_zts  *** 
     254      !!  
     255      !! **  Purpose :   Compute the now trend due to total advection of  
     256      !!       tracers and add it to the general trend of tracer equations 
     257      !! 
     258      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with 
     259      !!       corrected flux (monotonic correction). This version use sub- 
     260      !!       timestepping for the vertical advection which increases stability 
     261      !!       when vertical metrics are small. 
     262      !!       note: - this advection scheme needs a leap-frog time scheme 
     263      !! 
     264      !! ** Action : - update (pta) with the now advective tracer trends 
     265      !!             - save the trends  
     266      !!---------------------------------------------------------------------- 
     267      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace 
     268      ! 
     269      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
     270      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
     271      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     272      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
     273      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step 
     274      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     275      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
     276      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
     277      ! 
     278      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection 
     279      REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep 
     280      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices   
     281      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection 
     282      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps 
     283      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps 
     284      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection 
     285      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
     286      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
     287      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      - 
     288      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx_sav , zwy_sav 
     289      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 
     290      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 
     291      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 
     292      !!---------------------------------------------------------------------- 
     293      ! 
     294      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts') 
     295      ! 
     296      CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 
     297      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 
     298      CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 
     299      ! 
     300      IF( kt == kit000 )  THEN 
     301         IF(lwp) WRITE(numout,*) 
     302         IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 
     303         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     304      ENDIF 
     305      ! 
     306      l_trd = .FALSE. 
     307      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
     308      ! 
     309      IF( l_trd )  THEN 
     310         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     311         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp 
     312      ENDIF 
     313      ! 
     314      zwi(:,:,:) = 0._wp 
     315      z_rzts = 1._wp / REAL( jnzts, wp ) 
     316      zr_p2dt(:) = 1._wp / p2dt(:) 
     317      ! 
     318      !                                                          ! =========== 
     319      DO jn = 1, kjpt                                            ! tracer loop 
     320         !                                                       ! =========== 
     321         ! 1. Bottom value : flux set to zero 
     322         ! ---------------------------------- 
     323         zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
     324         zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
     325 
     326         ! 2. upstream advection with initial mass fluxes & intermediate update 
     327         ! -------------------------------------------------------------------- 
     328         ! upstream tracer flux in the i and j direction 
     329         DO jk = 1, jpkm1 
     330            DO jj = 1, jpjm1 
     331               DO ji = 1, fs_jpim1   ! vector opt. 
     332                  ! upstream scheme 
     333                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     334                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
     335                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
     336                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
     337                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) ) 
     338                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) ) 
     339               END DO 
     340            END DO 
     341         END DO 
     342 
     343         ! upstream tracer flux in the k direction 
     344         ! Surface value 
     345         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable 
     346         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
     347         ENDIF 
     348         ! Interior value 
     349         DO jk = 2, jpkm1 
     350            DO jj = 1, jpj 
     351               DO ji = 1, jpi 
     352                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     353                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     354                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 
     355               END DO 
     356            END DO 
     357         END DO 
     358 
     359         ! total advective trend 
     360         DO jk = 1, jpkm1 
     361            z2dtt = p2dt(jk) 
     362            DO jj = 2, jpjm1 
     363               DO ji = fs_2, fs_jpim1   ! vector opt. 
     364                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     365                  ! total intermediate advective trends 
     366                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     367                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     368                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     369                  ! update and guess with monotonic sheme 
     370                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
     371                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     372               END DO 
     373            END DO 
     374         END DO 
     375         !                             ! Lateral boundary conditions on zwi  (unchanged sign) 
     376         CALL lbc_lnk( zwi, 'T', 1. )   
     377 
     378         !                                 ! trend diagnostics (contribution of upstream fluxes) 
     379         IF( l_trd )  THEN  
     380            ! store intermediate advective trends 
     381            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
     382         END IF 
     383         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     384         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
     385           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
     386           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     387         ENDIF 
     388 
     389         ! 3. antidiffusive flux : high order minus low order 
     390         ! -------------------------------------------------- 
     391         ! antidiffusive flux on i and j 
     392 
     393 
     394         DO jk = 1, jpkm1 
     395 
     396            DO jj = 1, jpjm1 
     397               DO ji = 1, fs_jpim1   ! vector opt. 
     398                  zwx_sav(ji,jj) = zwx(ji,jj,jk) 
     399                  zwy_sav(ji,jj) = zwy(ji,jj,jk) 
     400 
     401                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 
     402                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 
     403               END DO 
     404            END DO 
     405 
     406            DO jj = 2, jpjm1         ! partial horizontal divergence 
     407               DO ji = fs_2, fs_jpim1 
     408                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
     409                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
     410               END DO 
     411            END DO 
     412 
     413            DO jj = 1, jpjm1 
     414               DO ji = 1, fs_jpim1   ! vector opt. 
     415                  zwx(ji,jj,jk) = zwx(ji,jj,jk)  - zwx_sav(ji,jj) 
     416                  zwy(ji,jj,jk) = zwy(ji,jj,jk)  - zwy_sav(ji,jj) 
     417               END DO 
     418            END DO 
     419         END DO 
     420       
     421         ! antidiffusive flux on k 
     422         zwz(:,:,1) = 0._wp        ! Surface value 
     423         zwz_sav(:,:,:) = zwz(:,:,:) 
     424         ! 
     425         ztrs(:,:,:,1) = ptb(:,:,:,jn) 
     426         zwzts(:,:,:) = 0._wp 
     427 
     428         DO jl = 1, jnzts                   ! Start of sub timestepping loop 
     429 
     430            IF( jl == 1 ) THEN              ! Euler forward to kick things off 
     431              jtb = 1   ;   jtn = 1   ;   jta = 2 
     432              zts(:) = p2dt(:) * z_rzts 
     433              jtaken = MOD( jnzts + 1 , 2)  ! Toggle to collect every second flux 
     434                                            ! starting at jl =1 if jnzts is odd;  
     435                                            ! starting at jl =2 otherwise 
     436            ELSEIF( jl == 2 ) THEN          ! First leapfrog step 
     437              jtb = 1   ;   jtn = 2   ;   jta = 3 
     438              zts(:) = 2._wp * p2dt(:) * z_rzts 
     439            ELSE                            ! Shuffle pointers for subsequent leapfrog steps 
     440              jtb = MOD(jtb,3) + 1 
     441              jtn = MOD(jtn,3) + 1 
     442              jta = MOD(jta,3) + 1 
     443            ENDIF 
     444            DO jk = 2, jpkm1          ! Interior value 
     445               DO jj = 2, jpjm1 
     446                  DO ji = fs_2, fs_jpim1 
     447                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 
     448                     IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux 
     449                  END DO 
     450               END DO 
     451            END DO 
     452 
     453            jtaken = MOD( jtaken + 1 , 2 ) 
     454 
     455            DO jk = 2, jpkm1          ! Interior value 
     456               DO jj = 2, jpjm1 
     457                  DO ji = fs_2, fs_jpim1 
     458                     zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     459                     ! total advective trends 
     460                     ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     461                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 
     462                  END DO 
     463               END DO 
     464            END DO 
     465 
     466         END DO 
     467 
     468         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping 
     469            DO jj = 2, jpjm1 
     470               DO ji = fs_2, fs_jpim1 
     471                  zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 
     472               END DO 
     473            END DO 
     474         END DO 
     475         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
     476         CALL lbc_lnk( zwz, 'W',  1. ) 
     477 
     478         ! 4. monotonicity algorithm 
     479         ! ------------------------- 
     480         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
     481 
     482 
     483         ! 5. final trend with corrected fluxes 
     484         ! ------------------------------------ 
     485         DO jk = 1, jpkm1 
     486            DO jj = 2, jpjm1 
     487               DO ji = fs_2, fs_jpim1   ! vector opt.   
     488                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     489                  ! total advective trends 
     490                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     491                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     492                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     493                  ! add them to the general tracer trends 
     494                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     495               END DO 
     496            END DO 
     497         END DO 
     498 
     499         !                                 ! trend diagnostics (contribution of upstream fluxes) 
     500         IF( l_trd )  THEN  
     501            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
     502            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     503            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
     504             
     505            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )    
     506            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
     507            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
     508         END IF 
     509         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     510         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
     511           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 
     512           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 
     513         ENDIF 
     514         ! 
     515      END DO 
     516      ! 
     517                   CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 
     518                   CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 
     519                   CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 
     520      IF( l_trd )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     521      ! 
     522      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts') 
     523      ! 
     524   END SUBROUTINE tra_adv_tvd_zts 
    249525 
    250526   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) 
Note: See TracChangeset for help on using the changeset viewer.