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 9019 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 – NEMO

Ignore:
Timestamp:
2017-12-13T15:58:53+01:00 (6 years ago)
Author:
timgraham
Message:

Merge of dev_CNRS_2017 into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r7753 r9019  
    99   !!---------------------------------------------------------------------- 
    1010   !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 
    11    !!  tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme  
    1211   !!                   with sub-time-stepping in the vertical direction 
    1312   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     
    2120   USE diaptr         ! poleward transport diagnostics 
    2221   USE diaar5         ! AR5 diagnostics 
    23    USE phycst, ONLY: rau0_rcp 
     22   USE phycst  , ONLY : rau0_rcp 
    2423   ! 
    2524   USE in_out_manager ! I/O manager 
    26    USE iom 
     25   USE iom            !  
    2726   USE lib_mpp        ! MPP library 
    2827   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    2928   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    30    USE wrk_nemo       ! Memory Allocation 
    3129   USE timing         ! Timing 
    3230 
     
    3432   PRIVATE 
    3533 
    36    PUBLIC   tra_adv_fct        ! routine called by traadv.F90 
    37    PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90 
    38    PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90 
     34   PUBLIC   tra_adv_fct        ! called by traadv.F90 
     35   PUBLIC   interp_4th_cpt     ! called by traadv_cen.F90 
    3936 
    4037   LOGICAL  ::   l_trd   ! flag to compute trends 
     
    5047#  include "vectopt_loop_substitute.h90" 
    5148   !!---------------------------------------------------------------------- 
    52    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     49   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    5350   !! $Id$ 
    5451   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7067      !! 
    7168      !! ** Action : - update pta  with the now advective tracer trends 
    72       !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T) 
     69      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T) 
    7370      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    7471      !!---------------------------------------------------------------------- 
     
    8885      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      - 
    8986      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      - 
    90       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
    91       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz, zptry 
    92       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d 
    93       !!---------------------------------------------------------------------- 
    94       ! 
    95       IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct') 
    96       ! 
    97       CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     87      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
     88      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry 
     89      !!---------------------------------------------------------------------- 
     90      ! 
     91      IF( ln_timing )   CALL timing_start('tra_adv_fct') 
    9892      ! 
    9993      IF( kt == kit000 )  THEN 
     
    10397      ENDIF 
    10498      ! 
    105       l_trd = .FALSE. 
     99      l_trd = .FALSE.            ! set local switches 
    106100      l_hst = .FALSE. 
    107101      l_ptr = .FALSE. 
    108       IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE. 
    109       IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE.  
    110       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    111          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
     102      IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
     103      IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
     104      IF(   cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     105         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
    112106      ! 
    113107      IF( l_trd .OR. l_hst )  THEN 
    114          CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     108         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) ) 
    115109         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp 
    116110      ENDIF 
    117111      ! 
    118112      IF( l_ptr ) THEN   
    119          CALL wrk_alloc( jpi, jpj, jpk, zptry ) 
     113         ALLOCATE( zptry(jpi,jpj,jpk) ) 
    120114         zptry(:,:,:) = 0._wp 
    121115      ENDIF 
     
    184178         END IF 
    185179         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    186          IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)  
     180         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:)  
    187181         ! 
    188182         !        !==  anti-diffusive flux : high order minus low order  ==! 
     
    308302         END DO 
    309303         ! 
    310          IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    311             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    312             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    313             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
     304         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
     305            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes  
     306            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes 
     307            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! 
     308            ! 
     309            IF( l_trd ) THEN              ! trend diagnostics 
     310               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
     311               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
     312               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
     313            ENDIF 
     314            !                             ! heat/salt transport 
     315            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
     316            ! 
     317            DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
    314318         ENDIF 
    315             ! 
    316          IF( l_trd ) THEN  
    317             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
    318             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
    319             CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
    320             ! 
    321          END IF 
    322          !                                !  heat/salt transport 
    323          IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
    324  
    325          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    326          IF( l_ptr ) THEN   
    327             zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     319         IF( l_ptr ) THEN              ! "Poleward" transports 
     320            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes 
    328321            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
     322            DEALLOCATE( zptry ) 
    329323         ENDIF 
    330324         ! 
    331325      END DO                     ! end of tracer loop 
    332326      ! 
    333                               CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
    334       IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    335       IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 
    336       ! 
    337       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct') 
     327      IF( ln_timing )   CALL timing_stop('tra_adv_fct') 
    338328      ! 
    339329   END SUBROUTINE tra_adv_fct 
    340  
    341  
    342    SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      & 
    343       &                                                  ptb, ptn, pta, kjpt, kn_fct_zts ) 
    344       !!---------------------------------------------------------------------- 
    345       !!                  ***  ROUTINE tra_adv_fct_zts  *** 
    346       !!  
    347       !! **  Purpose :   Compute the now trend due to total advection of  
    348       !!       tracers and add it to the general trend of tracer equations 
    349       !! 
    350       !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with 
    351       !!       corrected flux (monotonic correction). This version use sub- 
    352       !!       timestepping for the vertical advection which increases stability 
    353       !!       when vertical metrics are small. 
    354       !!       note: - this advection scheme needs a leap-frog time scheme 
    355       !! 
    356       !! ** Action : - update (pta) with the now advective tracer trends 
    357       !!             - save the trends  
    358       !!---------------------------------------------------------------------- 
    359       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    360       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    361       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    362       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    363       INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps 
    364       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    365       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    366       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
    367       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    368       ! 
    369       REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection 
    370       REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep 
    371       INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices   
    372       INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps 
    373       INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps 
    374       REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection 
    375       REAL(wp) ::   ztra            ! local scalar 
    376       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
    377       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      - 
    378       REAL(wp), POINTER, DIMENSION(:,:  )   ::   zwx_sav , zwy_sav 
    379       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav 
    380       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz 
    381       REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry 
    382       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs 
    383       !!---------------------------------------------------------------------- 
    384       ! 
    385       IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts') 
    386       ! 
    387       CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
    388       CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
    389       CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
    390       ! 
    391       IF( kt == kit000 )  THEN 
    392          IF(lwp) WRITE(numout,*) 
    393          IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype 
    394          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    395       ENDIF 
    396       ! 
    397       l_trd = .FALSE. 
    398       l_hst = .FALSE. 
    399       l_ptr = .FALSE. 
    400       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
    401       IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
    402       IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    403          &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
    404       ! 
    405       IF( l_trd .OR. l_hst )  THEN 
    406          CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    407          ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp 
    408       ENDIF 
    409       ! 
    410       IF( l_ptr ) THEN   
    411          CALL wrk_alloc( jpi, jpj,jpk, zptry ) 
    412          zptry(:,:,:) = 0._wp 
    413       ENDIF 
    414       zwi(:,:,:) = 0._wp 
    415       z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 
    416       zr_p2dt = 1._wp / p2dt 
    417       ! 
    418       ! surface & Bottom value : flux set to zero for all tracers 
    419       zwz(:,:, 1 ) = 0._wp 
    420       zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp 
    421       zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp 
    422       ! 
    423       !                                                          ! =========== 
    424       DO jn = 1, kjpt                                            ! tracer loop 
    425          !                                                       ! =========== 
    426          ! 
    427          ! Upstream advection with initial mass fluxes & intermediate update 
    428          DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction 
    429             DO jj = 1, jpjm1 
    430                DO ji = 1, fs_jpim1   ! vector opt. 
    431                   ! upstream scheme 
    432                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
    433                   zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    434                   zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    435                   zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    436                   zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) ) 
    437                   zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) ) 
    438                END DO 
    439             END DO 
    440          END DO 
    441          !                       ! upstream tracer flux in the k direction 
    442          DO jk = 2, jpkm1              ! Interior value 
    443             DO jj = 1, jpj 
    444                DO ji = 1, jpi 
    445                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    446                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    447                   zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    448                END DO 
    449             END DO 
    450          END DO 
    451          IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask) 
    452             IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value 
    453                DO jj = 1, jpj 
    454                   DO ji = 1, jpi 
    455                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)  
    456                   END DO 
    457                END DO    
    458             ELSE                             ! no cavities, surface value 
    459                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    460             ENDIF 
    461          ENDIF 
    462          ! 
    463          DO jk = 1, jpkm1         ! total advective trend 
    464             DO jj = 2, jpjm1 
    465                DO ji = fs_2, fs_jpim1   ! vector opt. 
    466                   !                             ! total intermediate advective trends 
    467                   ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    468                      &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    469                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj) 
    470                   !                             ! update and guess with monotonic sheme 
    471                   pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
    472                   zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    473                END DO 
    474             END DO 
    475          END DO 
    476          !                            
    477          CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign) 
    478          !                 
    479          IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes) 
    480             ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
    481          END IF 
    482          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    483          IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
    484  
    485          ! 3. anti-diffusive flux : high order minus low order 
    486          ! --------------------------------------------------- 
    487  
    488          DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes 
    489             ! 
    490             DO jj = 1, jpjm1 
    491                DO ji = 1, fs_jpim1   ! vector opt. 
    492                   zwx_sav(ji,jj) = zwx(ji,jj,jk) 
    493                   zwy_sav(ji,jj) = zwy(ji,jj,jk) 
    494                   ! 
    495                   zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 
    496                   zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 
    497                END DO 
    498             END DO 
    499             ! 
    500             DO jj = 2, jpjm1                    ! partial horizontal divergence 
    501                DO ji = fs_2, fs_jpim1 
    502                   zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    503                      &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    504                END DO 
    505             END DO 
    506             ! 
    507             DO jj = 1, jpjm1 
    508                DO ji = 1, fs_jpim1   ! vector opt. 
    509                   zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 
    510                   zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 
    511                END DO 
    512             END DO 
    513          END DO 
    514          ! 
    515          !                                !* vertical anti-diffusive flux 
    516          zwz_sav(:,:,:)   = zwz(:,:,:) 
    517          ztrs   (:,:,:,1) = ptb(:,:,:,jn) 
    518          ztrs   (:,:,1,2) = ptb(:,:,1,jn) 
    519          ztrs   (:,:,1,3) = ptb(:,:,1,jn) 
    520          zwzts  (:,:,:)   = 0._wp 
    521          ! 
    522          DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop 
    523             ! 
    524             IF( jl == 1 ) THEN                        ! Euler forward to kick things off 
    525                jtb = 1   ;   jtn = 1   ;   jta = 2 
    526                zts(:) = p2dt * z_rzts 
    527                jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux 
    528                !                                            ! starting at jl =1 if kn_fct_zts is odd;  
    529                !                                            ! starting at jl =2 otherwise 
    530             ELSEIF( jl == 2 ) THEN                    ! First leapfrog step 
    531                jtb = 1   ;   jtn = 2   ;   jta = 3 
    532                zts(:) = 2._wp * p2dt * z_rzts 
    533             ELSE                                      ! Shuffle pointers for subsequent leapfrog steps 
    534                jtb = MOD(jtb,3) + 1 
    535                jtn = MOD(jtn,3) + 1 
    536                jta = MOD(jta,3) + 1 
    537             ENDIF 
    538             DO jk = 2, jpkm1                          ! interior value 
    539                DO jj = 2, jpjm1 
    540                   DO ji = fs_2, fs_jpim1 
    541                      zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk) 
    542                      IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux 
    543                   END DO 
    544                END DO 
    545             END DO 
    546             IF( ln_linssh ) THEN                    ! top value (only in linear free surface case) 
    547                IF( ln_isfcav ) THEN                      ! ice-shelf cavities 
    548                   DO jj = 1, jpj 
    549                      DO ji = 1, jpi 
    550                         zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    551                      END DO 
    552                   END DO    
    553                ELSE                                      ! no ocean cavities 
    554                   zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
    555                ENDIF 
    556             ENDIF 
    557             ! 
    558             jtaken = MOD( jtaken + 1 , 2 ) 
    559             ! 
    560             DO jk = 2, jpkm1                             ! total advective trends 
    561                DO jj = 2, jpjm1 
    562                   DO ji = fs_2, fs_jpim1 
    563                      ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 & 
    564                         &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   & 
    565                         &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    566                   END DO 
    567                END DO 
    568             END DO 
    569             ! 
    570          END DO 
    571  
    572          DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping 
    573             DO jj = 2, jpjm1 
    574                DO ji = fs_2, fs_jpim1 
    575                   zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 
    576                END DO 
    577             END DO 
    578          END DO 
    579          CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    580          CALL lbc_lnk( zwz, 'W',  1. ) 
    581  
    582          ! 4. monotonicity algorithm 
    583          ! ------------------------- 
    584          CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
    585  
    586  
    587          ! 5. final trend with corrected fluxes 
    588          ! ------------------------------------ 
    589          DO jk = 1, jpkm1 
    590             DO jj = 2, jpjm1 
    591                DO ji = fs_2, fs_jpim1   ! vector opt.   
    592                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       & 
    593                      &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   & 
    594                      &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    595                END DO 
    596             END DO 
    597          END DO 
    598  
    599         ! 
    600          IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    601             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    602             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    603             ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    604          ENDIF 
    605             ! 
    606          IF( l_trd ) THEN  
    607             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
    608             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
    609             CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
    610             ! 
    611          END IF 
    612          !                                             ! heat/salt transport 
    613          IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
    614  
    615          !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    616          IF( l_ptr ) THEN   
    617             zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    618             CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
    619          ENDIF 
    620          ! 
    621       END DO 
    622       ! 
    623                               CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
    624                               CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
    625                               CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
    626       IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    627       IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 
    628       ! 
    629       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts') 
    630       ! 
    631    END SUBROUTINE tra_adv_fct_zts 
    632330 
    633331 
     
    653351      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars 
    654352      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    655       REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo 
    656       !!---------------------------------------------------------------------- 
    657       ! 
    658       IF( nn_timing == 1 )  CALL timing_start('nonosc') 
    659       ! 
    660       CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 
     353      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo 
     354      !!---------------------------------------------------------------------- 
     355      ! 
     356      IF( ln_timing )   CALL timing_start('nonosc') 
    661357      ! 
    662358      zbig  = 1.e+40_wp 
     
    734430      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign) 
    735431      ! 
    736       CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 
    737       ! 
    738       IF( nn_timing == 1 )  CALL timing_stop('nonosc') 
     432      IF( ln_timing )   CALL timing_stop('nonosc') 
    739433      ! 
    740434   END SUBROUTINE nonosc 
Note: See TracChangeset for help on using the changeset viewer.