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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90 @ 10806

Last change on this file since 10806 was 10806, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.
  • Property svn:keywords set to Id
File size: 10.7 KB
RevLine 
[5770]1MODULE traadv_cen
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen  ***
[6140]4   !! Ocean  tracers:   advective trend (2nd/4th order centered)
[5770]5   !!======================================================================
6   !! History :  3.7  ! 2014-05  (G. Madec)  original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
[6140]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
[5770]12   !!----------------------------------------------------------------------
[6140]13   USE dom_oce        ! ocean space and time domain
14   USE eosbn2         ! equation of state
15   USE traadv_fct     ! acces to routine interp_4th_cpt
16   USE trd_oce        ! trends: ocean variables
17   USE trdtra         ! trends manager: tracers
18   USE diaptr         ! poleward transport diagnostics
[7646]19   USE diaar5         ! AR5 diagnostics
[5770]20   !
[6140]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
[5770]25
26   IMPLICIT NONE
27   PRIVATE
28
[9019]29   PUBLIC   tra_adv_cen   ! called by traadv.F90
[5770]30   
31   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
32
[9019]33   LOGICAL ::   l_trd   ! flag to compute trends
34   LOGICAL ::   l_ptr   ! flag to compute poleward transport
35   LOGICAL ::   l_hst   ! flag to compute heat/salt transport
[7646]36
[5770]37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
[9598]40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10068]41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
[5770]43   !!----------------------------------------------------------------------
44CONTAINS
45
[10806]46   SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pw,     &
[10802]47      &                                               pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v ) 
[5770]48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_adv_cen  ***
50      !!                 
51      !! ** Purpose :   Compute the now trend due to the advection of tracers
52      !!      and add it to the general trend of passive tracer equations.
53      !!
54      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme
55      !!               using now fields (leap-frog scheme).
56      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
57      !!                = 4  ==>> 4th order    -        -       -      -
58      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical
59      !!                = 4  ==>> 4th order COMPACT  scheme     -      -
60      !!
[10802]61      !! ** Action : - update pt_rhs  with the now advective tracer trends
[6140]62      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
63      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[5770]64      !!----------------------------------------------------------------------
65      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
66      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
[10802]67      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms
[5770]68      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
69      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
70      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
71      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
[10806]72      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components
[10802]73      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! now tracer fields
74      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend
[5770]75      !
76      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
77      INTEGER  ::   ierr             ! local integer
78      REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
79      REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
[9019]80      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, ztu, ztv, ztw
[5770]81      !!----------------------------------------------------------------------
82      !
83      IF( kt == kit000 )  THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
86         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
87      ENDIF
[9019]88      !                          ! set local switches
[7646]89      l_trd = .FALSE.
90      l_hst = .FALSE.
91      l_ptr = .FALSE.
92      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE.
93      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
94      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
95         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
96      !
[6140]97      !                   
98      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
99      zwz(:,:,jpk) = 0._wp
[5770]100      !
101      DO jn = 1, kjpt            !==  loop over the tracers  ==!
102         !
[6140]103         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --!
[5770]104         !
[6140]105         CASE(  2  )                         !* 2nd order centered
[5770]106            DO jk = 1, jpkm1
107               DO jj = 1, jpjm1
108                  DO ji = 1, fs_jpim1   ! vector opt.
[10802]109                     zwx(ji,jj,jk) = 0.5_wp * pu(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji+1,jj  ,jk,jn) )
110                     zwy(ji,jj,jk) = 0.5_wp * pv(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji  ,jj+1,jk,jn) )
[5770]111                  END DO
112               END DO
113            END DO
114            !
[6140]115         CASE(  4  )                         !* 4th order centered
116            ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
[5770]117            ztv(:,:,jpk) = 0._wp
[6140]118            DO jk = 1, jpkm1                       ! masked gradient
119               DO jj = 2, jpjm1
[5770]120                  DO ji = fs_2, fs_jpim1   ! vector opt.
[10802]121                     ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
122                     ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
[5770]123                  END DO
124               END DO
125            END DO
[10425]126            CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond.
[5770]127            !
[6140]128            DO jk = 1, jpkm1                       ! Horizontal advective fluxes
[5770]129               DO jj = 2, jpjm1
130                  DO ji = 1, fs_jpim1   ! vector opt.
[10802]131                     zC2t_u = pt(ji,jj,jk,jn) + pt(ji+1,jj  ,jk,jn)   ! C2 interpolation of T at u- & v-points (x2)
132                     zC2t_v = pt(ji,jj,jk,jn) + pt(ji  ,jj+1,jk,jn)
[5770]133                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
134                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
135                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
136                     !                                                  ! C4 fluxes
[10802]137                     zwx(ji,jj,jk) =  0.5_wp * pu(ji,jj,jk) * zC4t_u
138                     zwy(ji,jj,jk) =  0.5_wp * pv(ji,jj,jk) * zC4t_v
[5770]139                  END DO
140               END DO
141            END DO         
142            !
143         CASE DEFAULT
144            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
145         END SELECT
146         !
[6140]147         SELECT CASE( kn_cen_v )       !--  Vertical fluxes  --!   (interior)
[5770]148         !
[6140]149         CASE(  2  )                         !* 2nd order centered
[5770]150            DO jk = 2, jpk
151               DO jj = 2, jpjm1
152                  DO ji = fs_2, fs_jpim1   ! vector opt.
[10806]153                     zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
[5770]154                  END DO
155               END DO
156            END DO
157            !
[6140]158         CASE(  4  )                         !* 4th order compact
[10802]159            CALL interp_4th_cpt( pt(:,:,:,jn) , ztw )      ! ztw = interpolated value of T at w-point
[5770]160            DO jk = 2, jpkm1
161               DO jj = 2, jpjm1
162                  DO ji = fs_2, fs_jpim1
[10806]163                     zwz(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
[5770]164                  END DO
165               END DO
166            END DO
167            !
168         END SELECT
169         !
[6140]170         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
[5770]171            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
172               DO jj = 1, jpj
173                  DO ji = 1, jpi
[10806]174                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
[5770]175                  END DO
176               END DO   
177            ELSE                                   ! no ice-shelf cavities (only ocean surface)
[10806]178               zwz(:,:,1) = pw(:,:,1) * pt(:,:,1,jn)
[5770]179            ENDIF
180         ENDIF
181         !               
182         DO jk = 1, jpkm1              !--  Divergence of advective fluxes  --!
183            DO jj = 2, jpjm1
184               DO ji = fs_2, fs_jpim1   ! vector opt.
[10802]185                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
[5770]186                     &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
187                     &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
[10802]188                     &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[5770]189               END DO
190            END DO
191         END DO
[6140]192         !                             ! trend diagnostics
[7646]193         IF( l_trd ) THEN
[10802]194            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu, pt(:,:,:,jn) )
195            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt(:,:,:,jn) )
[10806]196            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt(:,:,:,jn) )
[5770]197         END IF
[7646]198         !                                 ! "Poleward" heat and salt transports
[9019]199         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
[7646]200         !                                 !  heat and salt transport
[9019]201         IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
[5770]202         !
203      END DO
204      !
205   END SUBROUTINE tra_adv_cen
206   
207   !!======================================================================
208END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.