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

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

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

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.
File size: 11.0 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   USE diaar5         ! AR5 diagnostics
21   !
22   USE in_out_manager ! I/O manager
23   USE iom            ! IOM library
24   USE trc_oce        ! share passive tracers/Ocean variables
25   USE lib_mpp        ! MPP library
26   USE wrk_nemo       ! Memory Allocation
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_cen       ! routine called by step.F90
33   
34   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
35
36   LOGICAL :: l_trd   ! flag to compute trends
37   LOGICAL :: l_ptr   ! flag to compute poleward transport
38   LOGICAL :: l_hst   ! flag to compute heat/salt transport
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn,     &
50      &                                             ptn, pta, kjpt, 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 pta  with the now advective tracer trends
65      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
66      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
67      !!----------------------------------------------------------------------
68      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
69      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
70      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
71      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
72      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
73      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
74      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
75      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn             ! now tracer fields
76      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
77      !
78      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
79      INTEGER  ::   ierr             ! local integer
80      REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
81      REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
82      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zwy, zwz, ztu, ztv, ztw
83      !!----------------------------------------------------------------------
84      !
85      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen')
86      !
87      CALL wrk_alloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
88      !
89      IF( kt == kit000 )  THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
93      ENDIF
94      !
95      l_trd = .FALSE.
96      l_hst = .FALSE.
97      l_ptr = .FALSE.
98      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE.
99      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
100      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
101         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
102      !
103      !                   
104      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
105      zwz(:,:,jpk) = 0._wp
106      !
107      DO jn = 1, kjpt            !==  loop over the tracers  ==!
108         !
109         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --!
110         !
111         CASE(  2  )                         !* 2nd order centered
112            DO jk = 1, jpkm1
113               DO jj = 1, jpjm1
114                  DO ji = 1, fs_jpim1   ! vector opt.
115                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
116                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
117                  END DO
118               END DO
119            END DO
120            !
121         CASE(  4  )                         !* 4th order centered
122            ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
123            ztv(:,:,jpk) = 0._wp
124            DO jk = 1, jpkm1                       ! masked gradient
125               DO jj = 2, jpjm1
126                  DO ji = fs_2, fs_jpim1   ! vector opt.
127                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
128                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
129                  END DO
130               END DO
131            END DO
132            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
133            !
134            DO jk = 1, jpkm1                       ! Horizontal advective fluxes
135               DO jj = 2, jpjm1
136                  DO ji = 1, fs_jpim1   ! vector opt.
137                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
138                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
139                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
140                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
141                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
142                     !                                                  ! C4 fluxes
143                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u
144                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v
145                  END DO
146               END DO
147            END DO         
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 jk = 2, jpk
157               DO jj = 2, jpjm1
158                  DO ji = fs_2, fs_jpim1   ! vector opt.
159                     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)
160                  END DO
161               END DO
162            END DO
163            !
164         CASE(  4  )                         !* 4th order compact
165            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point
166            DO jk = 2, jpkm1
167               DO jj = 2, jpjm1
168                  DO ji = fs_2, fs_jpim1
169                     zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
170                  END DO
171               END DO
172            END DO
173            !
174         END SELECT
175         !
176         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
177            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
178               DO jj = 1, jpj
179                  DO ji = 1, jpi
180                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 
181                  END DO
182               END DO   
183            ELSE                                   ! no ice-shelf cavities (only ocean surface)
184               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
185            ENDIF
186         ENDIF
187         !               
188         DO jk = 1, jpkm1              !--  Divergence of advective fluxes  --!
189            DO jj = 2, jpjm1
190               DO ji = fs_2, fs_jpim1   ! vector opt.
191                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)    &
192                     &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
193                     &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
194                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
195               END DO
196            END DO
197         END DO
198         !                             ! trend diagnostics
199         IF( l_trd ) THEN
200            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
201            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
202            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
203         END IF
204         !                                 ! "Poleward" heat and salt transports
205         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
206         !                                 !  heat and salt transport
207         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
208         !
209      END DO
210      !
211      CALL wrk_dealloc( jpi,jpj,jpk,   zwx, zwy, zwz, ztu, ztv, ztw )
212      !
213      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen')
214      !
215   END SUBROUTINE tra_adv_cen
216   
217   !!======================================================================
218END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.