source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90 @ 5836

Last change on this file since 5836 was 5836, checked in by cetlod, 5 years ago

merge the simplification branch onto the trunk, see ticket #1612

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