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.
traadv_cen.F90 in branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 8 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

File size: 10.7 KB
Line 
1MODULE traadv_cen
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order centered)
5   !!======================================================================
6   !! History :  3.7  ! 2014-05  (G. Madec)  original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_cen : update the tracer trend with the advection trends using a centered or scheme (2nd or 4th order)
11   !!                 NB: on the vertical it is actually a 4th order COMPACT scheme which is used
12   !!----------------------------------------------------------------------
13   USE oce, ONLY: tsn  ! now ocean temperature and salinity
14   USE dom_oce         ! ocean space and time domain
15   USE eosbn2          ! equation of state
16   USE traadv_fct      ! acces to routine interp_4th_cpt
17   USE trd_oce         ! trends: ocean variables
18   USE trdtra          ! trends manager: tracers
19   USE diaptr          ! poleward transport diagnostics
20   !
21   USE in_out_manager  ! I/O manager
22   USE iom             ! IOM library
23   USE trc_oce         ! share passive tracers/Ocean variables
24   USE lib_mpp         ! MPP library
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   tra_adv_cen       ! routine called by step.F90
32   
33   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
39   !! $Id: traadv_cen2.F90 5737 2015-09-13 07:42:41Z gm $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn,     &
45      &                                             ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tra_adv_cen  ***
48      !!                 
49      !! ** Purpose :   Compute the now trend due to the advection of tracers
50      !!      and add it to the general trend of passive tracer equations.
51      !!
52      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme
53      !!               using now fields (leap-frog scheme).
54      !!
55      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
56      !!                = 4  ==>> 4th order    -        -       -      -
57      !!
58      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical
59      !!                = 4  ==>> 4th order COMPACT  scheme     -      -
60      !!
61      !! ** Action :  - update pta  with the now advective tracer trends
62      !!              - send trends to trdtra module for further diagnostcs
63      !!----------------------------------------------------------------------
64      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
65      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
66      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
67      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
68      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
69      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
70      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
71      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! now tracer fields
72      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
73      !
74      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
75      INTEGER  ::   ierr             ! local integer
76      REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
77      REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zwy, zwz, ztu, ztv, ztw
79      !!----------------------------------------------------------------------
80      !
81      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen')
82      !
83      CALL wrk_alloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
84      !
85      IF( kt == kit000 )  THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
89      ENDIF
90      !
91      !                          ! surface & bottom values
92      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all
93                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface
94      !
95      DO jn = 1, kjpt            !==  loop over the tracers  ==!
96         !
97         SELECT CASE( kn_cen_h )          !--  Horizontal fluxes  --!
98         !
99         CASE(  2  )                               ! 2nd order centered
100            DO jk = 1, jpkm1
101               DO jj = 1, jpjm1
102                  DO ji = 1, fs_jpim1   ! vector opt.
103                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
104                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
105                  END DO
106               END DO
107            END DO
108            !
109         CASE(  4  )                               ! 4th order centered
110            ztu(:,:,jpk) = 0._wp                            ! Bottom value : flux set to zero
111            ztv(:,:,jpk) = 0._wp
112            DO jk = 1, jpkm1                                 ! gradient
113               DO jj = 2, jpjm1                                   ! masked derivative
114                  DO ji = fs_2, fs_jpim1   ! vector opt.
115                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
116                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
117                  END DO
118               END DO
119            END DO
120            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
121            !
122            DO jk = 1, jpkm1                                 ! Horizontal advective fluxes
123               DO jj = 2, jpjm1
124                  DO ji = 1, fs_jpim1   ! vector opt.
125                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
126                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
127                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
128                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
129                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
130                     !                                                  ! C4 fluxes
131                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u
132                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v
133                  END DO
134               END DO
135            END DO         
136            !
137         CASE DEFAULT
138            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
139         END SELECT
140         !
141         !                             !==  Vertical fluxes  ==!
142         !
143         SELECT CASE( kn_cen_v )             !* interior fluxes
144         !
145         CASE(  2  )                               ! 2nd order centered
146            DO jk = 2, jpk
147               DO jj = 2, jpjm1
148                  DO ji = fs_2, fs_jpim1   ! vector opt.
149                     zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
150                  END DO
151               END DO
152            END DO
153            !
154         CASE(  4  )                               ! 4th order centered
155            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
156            DO jk = 2, jpkm1
157               DO jj = 2, jpjm1
158                  DO ji = fs_2, fs_jpim1
159                     zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
160                  END DO
161               END DO
162            END DO
163            !
164         END SELECT
165         !
166         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask)
167            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
168               DO jj = 1, jpj
169                  DO ji = 1, jpi
170                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
171                  END DO
172               END DO   
173            ELSE                                   ! no ice-shelf cavities (only ocean surface)
174               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
175            ENDIF
176         ENDIF
177         !               
178         DO jk = 1, jpkm1              !--  Divergence of advective fluxes  --!
179            DO jj = 2, jpjm1
180               DO ji = fs_2, fs_jpim1   ! vector opt.
181                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
182                     &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
183                     &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
184                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) / ( e1e2t(ji,jj) *  e3t_n(ji,jj,jk) )
185               END DO
186            END DO
187         END DO
188         !                                 ! trend diagnostics
189         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
190            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
191            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
192            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
193         END IF
194         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
195         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
196           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
197           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
198         ENDIF
199         !
200      END DO
201      !
202      CALL wrk_dealloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
203      !
204      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen')
205      !
206   END SUBROUTINE tra_adv_cen
207   
208   !!======================================================================
209END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.