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/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90 @ 7236

Last change on this file since 7236 was 7236, checked in by timgraham, 7 years ago

All changes related to diaptr (basin heat transports and transport components)
File size: 10.4 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 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      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
55      !!                = 4  ==>> 4th order    -        -       -      -
56      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical
57      !!                = 4  ==>> 4th order COMPACT  scheme     -      -
58      !!
59      !! ** Action : - update pta  with the now advective tracer trends
60      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
61      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
62      !!----------------------------------------------------------------------
63      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
64      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
65      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
66      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
67      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
68      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
69      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
70      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! now tracer fields
71      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
72      !
73      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
74      INTEGER  ::   ierr             ! local integer
75      REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
76      REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
77      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zwy, zwz, ztu, ztv, ztw
78      !!----------------------------------------------------------------------
79      !
80      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen')
81      !
82      CALL wrk_alloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
83      !
84      IF( kt == kit000 )  THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
87         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
88      ENDIF
89      !
90      !                   
91      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
92      zwz(:,:,jpk) = 0._wp
93      !
94      DO jn = 1, kjpt            !==  loop over the tracers  ==!
95         !
96         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --!
97         !
98         CASE(  2  )                         !* 2nd order centered
99            DO jk = 1, jpkm1
100               DO jj = 1, jpjm1
101                  DO ji = 1, fs_jpim1   ! vector opt.
102                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
103                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
104                  END DO
105               END DO
106            END DO
107            !
108         CASE(  4  )                         !* 4th order centered
109            ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
110            ztv(:,:,jpk) = 0._wp
111            DO jk = 1, jpkm1                       ! masked gradient
112               DO jj = 2, jpjm1
113                  DO ji = fs_2, fs_jpim1   ! vector opt.
114                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
115                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
116                  END DO
117               END DO
118            END DO
119            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
120            !
121            DO jk = 1, jpkm1                       ! Horizontal advective fluxes
122               DO jj = 2, jpjm1
123                  DO ji = 1, fs_jpim1   ! vector opt.
124                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
125                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
126                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
127                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
128                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
129                     !                                                  ! C4 fluxes
130                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u
131                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v
132                  END DO
133               END DO
134            END DO         
135            !
136         CASE DEFAULT
137            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
138         END SELECT
139         !
140         SELECT CASE( kn_cen_v )       !--  Vertical fluxes  --!   (interior)
141         !
142         CASE(  2  )                         !* 2nd order centered
143            DO jk = 2, jpk
144               DO jj = 2, jpjm1
145                  DO ji = fs_2, fs_jpim1   ! vector opt.
146                     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)
147                  END DO
148               END DO
149            END DO
150            !
151         CASE(  4  )                         !* 4th order compact
152            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point
153            DO jk = 2, jpkm1
154               DO jj = 2, jpjm1
155                  DO ji = fs_2, fs_jpim1
156                     zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
157                  END DO
158               END DO
159            END DO
160            !
161         END SELECT
162         !
163         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
164            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
165               DO jj = 1, jpj
166                  DO ji = 1, jpi
167                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 
168                  END DO
169               END DO   
170            ELSE                                   ! no ice-shelf cavities (only ocean surface)
171               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
172            ENDIF
173         ENDIF
174         !               
175         DO jk = 1, jpkm1              !--  Divergence of advective fluxes  --!
176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.
178                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
179                     &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
180                     &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
181                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
182               END DO
183            END DO
184         END DO
185         !                             ! trend diagnostics
186         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
187            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
188            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
189            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
190         END IF
191         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
192         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
193         !
194      END DO
195      !
196      CALL wrk_dealloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
197      !
198      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen')
199      !
200   END SUBROUTINE tra_adv_cen
201   
202   !!======================================================================
203END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.