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 NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_cen.F90 @ 13516

Last change on this file since 13516 was 13516, checked in by hadcv, 4 years ago

Tiling for tra_adv
  • Property svn:keywords set to Id
File size: 11.5 KB
Line 
1MODULE traadv_cen
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen  ***
4   !! Ocean  tracers:   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 dom_oce        ! ocean space and time domain
14   ! TEMP: This change not necessary after trd_tra is tiled
15   USE domain, ONLY : dom_tile
16   USE eosbn2         ! equation of state
17   USE traadv_fct     ! acces to routine interp_4th_cpt
18   USE trd_oce        ! trends: ocean variables
19   USE trdtra         ! trends manager: tracers
20   USE diaptr         ! poleward transport diagnostics
21   USE diaar5         ! AR5 diagnostics
22   !
23   USE in_out_manager ! I/O manager
24   USE iom            ! IOM library
25   USE trc_oce        ! share passive tracers/Ocean variables
26   USE lib_mpp        ! MPP library
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   tra_adv_cen   ! called by traadv.F90
32   
33   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
34
35   LOGICAL ::   l_trd   ! flag to compute trends
36   LOGICAL ::   l_ptr   ! flag to compute poleward transport
37   LOGICAL ::   l_hst   ! flag to compute heat/salt transport
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     &
50      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_cen  ***
53      !!                 
54      !! ** Purpose :   Compute the now trend due to the advection of tracers
55      !!      and add it to the general trend of passive tracer equations.
56      !!
57      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme
58      !!               using now fields (leap-frog scheme).
59      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
60      !!                = 4  ==>> 4th order    -        -       -      -
61      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical
62      !!                = 4  ==>> 4th order COMPACT  scheme     -      -
63      !!
64      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
65      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
66      !!             - poleward advective heat and salt transport (l_diaptr=T)
67      !!----------------------------------------------------------------------
68      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
69      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs       ! ocean time level indices
70      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
71      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
72      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
73      INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
74      INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
75      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
76      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
77      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
78      !
79      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
80      INTEGER  ::   ierr             ! local integer
81      ! TEMP: This change not necessary after trd_tra is tiled
82      INTEGER  ::   itile
83      REAL(wp) ::   zC2t_u, zC2t_v   ! local scalars
84      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   zwx, zwy, zwz, ztu, ztv, ztw, zltu, zltv
85      ! TEMP: This change not necessary after trd_tra is tiled
86      REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE ::   ztrdx, ztrdy, ztrdz
87      !!----------------------------------------------------------------------
88      ! TEMP: This change not necessary after trd_tra is tiled
89      itile = ntile
90      !
91      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
92         IF( kt == kit000 )  THEN
93            IF(lwp) WRITE(numout,*)
94            IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
95            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
96         ENDIF
97         !                          ! set local switches
98         l_trd = .FALSE.
99         l_hst = .FALSE.
100         l_ptr = .FALSE.
101         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
102         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE.
103         IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
104            &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
105
106         ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
107         IF( kt == kit000 .AND. l_trd ) THEN
108            ALLOCATE( ztrdx(jpi,jpj,jpk,jpts), ztrdy(jpi,jpj,jpk,jpts), ztrdz(jpi,jpj,jpk,jpts) )
109         ENDIF
110      ENDIF
111      !
112      !                   
113      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
114      zwz(:,:,jpk) = 0._wp
115      !
116      DO jn = 1, kjpt            !==  loop over the tracers  ==!
117         !
118         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --!
119         !
120         CASE(  2  )                         !* 2nd order centered
121            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
122               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
123               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
124            END_3D
125            !
126         CASE(  4  )                         !* 4th order centered
127            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
128            zltv(:,:,jpk) = 0._wp
129            DO jk = 1, jpkm1                 ! Laplacian
130               DO_2D( 1, 0, 1, 0 )
131                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
132                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
133               END_2D
134               DO_2D( 0, 0, 0, 0 )
135                  zltu(ji,jj,jk) = ztu(ji,jj,jk) + ztu(ji-1,jj,jk)
136                  zltv(ji,jj,jk) = ztv(ji,jj,jk) + ztv(ji,jj-1,jk)
137               END_2D
138            END DO
139            CALL lbc_lnk_multi( 'traadv_cen', zltu, 'T', 1. , zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
140            !
141            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
142               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
143               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
144               !                                                  ! C4 fluxes
145               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + r1_6 * (zltu(ji,jj,jk) - zltu(ji+1,jj,jk)) )
146               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + r1_6 * (zltv(ji,jj,jk) - zltv(ji,jj+1,jk)) )
147            END_3D
148            !
149         CASE DEFAULT
150            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
151         END SELECT
152         !
153         SELECT CASE( kn_cen_v )       !--  Vertical fluxes  --!   (interior)
154         !
155         CASE(  2  )                         !* 2nd order centered
156            DO_3D( 0, 0, 0, 0, 2, jpk )
157               zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
158            END_3D
159            !
160         CASE(  4  )                         !* 4th order compact
161            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )      ! ztw = interpolated value of T at w-point
162            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
163               zwz(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
164            END_3D
165            !
166         END SELECT
167         !
168         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
169            ! TODO: NOT TESTED- requires isf
170            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
171               DO_2D( 1, 1, 1, 1 )
172                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
173               END_2D
174            ELSE                                   ! no ice-shelf cavities (only ocean surface)
175               DO_2D( 1, 1, 1, 1 )
176                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
177               END_2D
178            ENDIF
179         ENDIF
180         !               
181         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
182            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
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)  ) &
186               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
187         END_3D
188         !                             ! trend diagnostics
189         ! TEMP: These changes not necessary after trd_tra is tiled
190         IF( l_trd ) THEN
191            DO_3D( 1, 0, 1, 0, 1, jpk )
192               ztrdx(ji,jj,jk,jn) = zwx(ji,jj,jk)
193               ztrdy(ji,jj,jk,jn) = zwy(ji,jj,jk)
194               ztrdz(ji,jj,jk,jn) = zwz(ji,jj,jk)
195            END_3D
196
197            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
198               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
199
200               ! TODO: TO BE TILED- trd_tra
201               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx(:,:,:,jn), pU, pt(:,:,:,jn,Kmm) )
202               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy(:,:,:,jn), pV, pt(:,:,:,jn,Kmm) )
203               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz(:,:,:,jn), pW, pt(:,:,:,jn,Kmm) )
204
205               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
206            ENDIF
207         ENDIF
208         !                                 ! "Poleward" heat and salt transports
209         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
210         !                                 !  heat and salt transport
211         IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
212         !
213      END DO
214      !
215   END SUBROUTINE tra_adv_cen
216   
217   !!======================================================================
218END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.