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_mus.F90 in NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/TRA/traadv_mus.F90 @ 11719

Last change on this file since 11719 was 11719, checked in by francesca, 5 years ago

add extra halo support- ticket #2009
  • Property svn:keywords set to Id
File size: 21.5 KB
Line 
1MODULE traadv_mus
2   !!======================================================================
3   !!                       ***  MODULE  traadv_mus  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :       !  2000-06  (A.Estublier)  for passive tracers
7   !!                 !  2001-08  (E.Durand, G.Madec)  adapted for T & S
8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
10   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
11   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_adv_mus   : update the tracer trend with the horizontal
16   !!                   and vertical advection trends using MUSCL scheme
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE trc_oce        ! share passive tracers/Ocean variables
20   USE dom_oce        ! ocean space and time domain
21   USE trd_oce        ! trends: ocean variables
22   USE trdtra         ! tracers trends manager
23   USE sbcrnf         ! river runoffs
24   USE diaptr         ! poleward transport diagnostics
25   USE diaar5         ! AR5 diagnostics
26
27   !
28   USE iom            ! XIOS library
29   USE in_out_manager ! I/O manager
30   USE lib_mpp        ! distribued memory computing
31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
33   USE halo_mng
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_adv_mus   ! routine called by traadv.F90
39   
40   REAL(wp), ALLOCATABLE, DIMENSION(:,:    ) ::   r1_e1e2t_exh2, r1_e1e2u_exh2, r1_e1e2v_exh2
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:    ) ::   rnfmsk_exh2, upsmsk_exh2, mikt_exh2
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) ::   tmask_exh2, wmask_exh2, umask_exh2, vmask_exh2
43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) ::   e3u_n_exh2, e3v_n_exh2, e3t_n_exh2, e3w_n_exh2
44   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) ::   pun_exh2, pvn_exh2, pwn_exh2   ! 3 ocean velocity components
45   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ptb_exh2, pta_exh2        ! before and now tracer fields
46
47   
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
49   !                                                           !  and in closed seas (orca 2 and 1 configurations)
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
51   
52   LOGICAL  ::   l_trd   ! flag to compute trends
53   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
54   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
55   
56   INTEGER :: jphls = 2
57
58   !! * Substitutions
59#  include "vectopt_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id$
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             &
68      &                                              ptb, pta, kjpt, ld_msc_ups )
69      !!----------------------------------------------------------------------
70      !!                    ***  ROUTINE tra_adv_mus  ***
71      !!
72      !! ** Purpose :   Compute the now trend due to total advection of tracers
73      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
74      !!              Conservation Laws) and add it to the general tracer trend.
75      !!
76      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
77      !!              ld_msc_ups=T :
78      !!
79      !! ** Action : - update pta  with the now advective tracer trends
80      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
81      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
82      !!
83      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
84      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
85      !!----------------------------------------------------------------------
86      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
87      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
90      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
91      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
95      !
96      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
97      INTEGER  ::   last_khls, ierr             ! local integer
98      REAL(wp) ::   zu, z0u, zzwx, zw , zalpha   ! local scalars
99      REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
100      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwx, zslpx   ! 3D workspace
101      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwy, zslpy   ! -      -
102      !!----------------------------------------------------------------------
103      !
104     
105   CALL halo_mng_set(jphls)
106   
107   ALLOCATE(zwx(jplbi:jpi,jplbj:jpj,jpk))
108   ALLOCATE(zwy(jplbi:jpi,jplbj:jpj,jpk))
109   ALLOCATE(zslpx(jplbi:jpi,jplbj:jpj,jpk))
110   ALLOCATE(zslpy(jplbi:jpi,jplbj:jpj,jpk))
111
112   IF (kt==kit000) THEN
113      if (.not. allocated(pun_exh2))   ALLOCATE(pun_exh2(jplbi:jpi,jplbj:jpj,jpk))
114      if (.not. allocated(pvn_exh2))   ALLOCATE(pvn_exh2(jplbi:jpi,jplbj:jpj,jpk))
115      if (.not. allocated(pwn_exh2))   ALLOCATE(pwn_exh2(jplbi:jpi,jplbj:jpj,jpk))
116      if (.not. allocated(ptb_exh2))   ALLOCATE(ptb_exh2(jplbi:jpi,jplbj:jpj,jpk,kjpt))
117      if (.not. allocated(pta_exh2))   ALLOCATE(pta_exh2(jplbi:jpi,jplbj:jpj,jpk,kjpt))
118      if (.not. allocated(r1_e1e2t_exh2)) ALLOCATE(r1_e1e2t_exh2(jplbi:jpi,jplbj:jpj))
119      if (.not. allocated(r1_e1e2u_exh2)) ALLOCATE(r1_e1e2u_exh2(jplbi:jpi,jplbj:jpj))
120      if (.not. allocated(r1_e1e2v_exh2)) ALLOCATE(r1_e1e2v_exh2(jplbi:jpi,jplbj:jpj))
121      if (.not. allocated(tmask_exh2)) ALLOCATE(tmask_exh2(jplbi:jpi,jplbj:jpj,jpk))
122      if (.not. allocated(wmask_exh2)) ALLOCATE(wmask_exh2(jplbi:jpi,jplbj:jpj,jpk))
123      if (.not. allocated(umask_exh2)) ALLOCATE(umask_exh2(jplbi:jpi,jplbj:jpj,jpk))
124      if (.not. allocated(vmask_exh2)) ALLOCATE(vmask_exh2(jplbi:jpi,jplbj:jpj,jpk))
125      if (.not. allocated(e3u_n_exh2)) ALLOCATE(e3u_n_exh2(jplbi:jpi,jplbj:jpj,jpk))
126      if (.not. allocated(e3v_n_exh2)) ALLOCATE(e3v_n_exh2(jplbi:jpi,jplbj:jpj,jpk))
127      if (.not. allocated(e3t_n_exh2)) ALLOCATE(e3t_n_exh2(jplbi:jpi,jplbj:jpj,jpk))
128      if (.not. allocated(e3w_n_exh2)) ALLOCATE(e3w_n_exh2(jplbi:jpi,jplbj:jpj,jpk))
129      IF( ln_isfcav.and..not.allocated(mikt_exh2)) ALLOCATE(mikt_exh2(jplbi:jpi,jplbj:jpj))
130      IF( ld_msc_ups.and..not.allocated(rnfmsk_exh2)) ALLOCATE(rnfmsk_exh2(jplbi:jpi,jplbj:jpj))
131      IF( ld_msc_ups.and..not.allocated(upsmsk_exh2)) ALLOCATE(upsmsk_exh2(jplbi:jpi,jplbj:jpj))
132
133      CALL halo_mng_copy(r1_e1e2t, r1_e1e2t_exh2)
134      CALL halo_mng_copy(r1_e1e2u, r1_e1e2u_exh2)
135      CALL halo_mng_copy(r1_e1e2v, r1_e1e2v_exh2)
136      CALL halo_mng_copy(tmask, tmask_exh2)
137      CALL halo_mng_copy(wmask, wmask_exh2)
138      CALL halo_mng_copy(umask, umask_exh2)
139      CALL halo_mng_copy(vmask, vmask_exh2)
140
141      CALL lbc_lnk( 'traadv_mus', r1_e1e2u_exh2, 'U', -1.)
142      CALL lbc_lnk( 'traadv_mus', r1_e1e2v_exh2, 'V', -1.)
143      CALL lbc_lnk( 'traadv_mus', r1_e1e2t_exh2, 'T', 1.)
144      CALL lbc_lnk( 'traadv_mus', tmask_exh2, 'T', 1. )
145      CALL lbc_lnk( 'traadv_mus', wmask_exh2, 'W', 1.)
146      CALL lbc_lnk( 'traadv_mus', umask_exh2, 'U', 1. )
147      CALL lbc_lnk( 'traadv_mus', vmask_exh2, 'V', 1.)
148   ENDIF
149   
150   IF( ln_isfcav ) THEN ; CALL halo_mng_copy(REAL(mikt), mikt_exh2) ; CALL lbc_lnk( 'traadv_mus', mikt_exh2, 'T', 1.) ; ENDIF
151   IF( ld_msc_ups) THEN ; CALL halo_mng_copy(rnfmsk, rnfmsk_exh2) ; CALL lbc_lnk( 'traadv_mus', rnfmsk_exh2, 'T', 1.) ; ENDIF
152   IF( ld_msc_ups) THEN ; CALL halo_mng_copy(upsmsk, upsmsk_exh2) ; CALL lbc_lnk( 'traadv_mus', upsmsk_exh2, 'T', 1.) ; ENDIF
153
154   CALL halo_mng_copy(e3u_n, e3u_n_exh2)
155   CALL halo_mng_copy(e3v_n, e3v_n_exh2)
156   CALL halo_mng_copy(e3t_n, e3t_n_exh2)
157   CALL halo_mng_copy(e3w_n, e3w_n_exh2)
158   CALL halo_mng_copy(pun, pun_exh2)
159   CALL halo_mng_copy(pvn, pvn_exh2)
160   CALL halo_mng_copy(pwn, pwn_exh2)
161   CALL halo_mng_copy(ptb, ptb_exh2)
162   CALL halo_mng_copy(pta, pta_exh2)
163
164   CALL lbc_lnk( 'traadv_mus', e3u_n_exh2, 'U', -1., pfillval = 1.0_wp )
165   CALL lbc_lnk( 'traadv_mus', e3v_n_exh2, 'V', -1., pfillval = 1.0_wp )
166   CALL lbc_lnk( 'traadv_mus', e3t_n_exh2, 'T', 1., pfillval = 1.0_wp )
167   CALL lbc_lnk( 'traadv_mus', e3w_n_exh2, 'W', 1., pfillval = 1.0_wp )
168   CALL lbc_lnk( 'traadv_mus', pun_exh2, 'U', -1.)
169   CALL lbc_lnk( 'traadv_mus', pvn_exh2, 'V', -1.)
170   CALL lbc_lnk( 'traadv_mus', pwn_exh2, 'W', 1.)
171   CALL lbc_lnk( 'traadv_mus', pta_exh2, 'T', 1.)
172   CALL lbc_lnk( 'traadv_mus', ptb_exh2, 'T', 1.)
173
174#     define pun pun_exh2
175#     define pvn pvn_exh2
176#     define pwn pwn_exh2
177#     define ptb ptb_exh2
178#     define pta pta_exh2
179#     define r1_e1e2t r1_e1e2t_exh2
180#     define r1_e1e2u r1_e1e2u_exh2
181#     define r1_e1e2v r1_e1e2v_exh2
182#     define tmask tmask_exh2
183#     define wmask wmask_exh2
184#     define umask umask_exh2
185#     define vmask vmask_exh2
186#     define e3u_n e3u_n_exh2
187#     define e3v_n e3v_n_exh2
188#     define e3t_n e3t_n_exh2
189#     define e3w_n e3w_n_exh2
190#     define mikt mikt_exh2
191#     define rnfmsk rnfmsk_exh2
192#     define upsmsk upsmsk_exh2
193
194      IF( kt == kit000 )  THEN
195         IF(lwp) WRITE(numout,*)
196         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
197         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
198         IF(lwp) WRITE(numout,*) '~~~~~~~'
199         IF(lwp) WRITE(numout,*)
200         !
201         ! Upstream / MUSCL scheme indicator
202         !
203         ALLOCATE( xind(jplbi:jpi,jplbj:jpj,jpk), STAT=ierr )
204         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
205         !
206         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
207            ALLOCATE( upsmsk(jplbi:jpi,jplbj:jpj), STAT=ierr )
208            upsmsk(:,:) = 0._wp                             ! not upstream by default
209            !
210            DO jk = 1, jpkm1
211               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
212                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
213                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
214            END DO
215         ENDIF 
216         !
217      ENDIF 
218      !     
219      l_trd = .FALSE.
220      l_hst = .FALSE.
221      l_ptr = .FALSE.
222      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
223      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
224      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
225         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
226      !
227      DO jn = 1, kjpt            !==  loop over the tracers  ==!
228         !
229         !                          !* Horizontal advective fluxes
230         !
231         !                                !-- first guess of the slopes
232         zwx(:,:,jpk) = 0._wp                   ! bottom values
233         zwy(:,:,jpk) = 0._wp 
234         DO jk = 1, jpkm1                       ! interior values
235            DO jj = jplbj, jpj-1     
236               DO ji = jplbi, jpi-1   ! vector opt.
237                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
238                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
239               END DO
240           END DO
241         END DO
242         ! lateral boundary conditions   (changed sign)
243         !CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
244         !                                !-- Slopes of tracer
245         zslpx(:,:,jpk) = 0._wp                 ! bottom values
246         zslpy(:,:,jpk) = 0._wp
247         DO jk = 1, jpkm1                       ! interior values
248            DO jj = jplbj+1, jpj
249               DO ji = jplbi+1, jpi   ! vector opt.
250                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
251                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
252                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
253                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
254               END DO
255            END DO
256         END DO
257         !
258         DO jk = 1, jpkm1                 !-- Slopes limitation
259            DO jj = jplbj+1, jpj
260               DO ji = jplbi+1, jpi   ! vector opt.
261                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
262                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
263                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
264                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
265                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
266                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
267               END DO
268           END DO
269         END DO
270         !
271         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes
272            DO jj = jplbj+1, jpj-1
273               DO ji = jplbi+1, jpi-1   ! vector opt.
274                  ! MUSCL fluxes
275                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
276                  zalpha = 0.5 - z0u
277                  zu  = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
278                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
279                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
280                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
281                  !
282                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
283                  zalpha = 0.5 - z0v
284                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
285                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
286                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
287                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
288               END DO
289            END DO
290         END DO
291         !CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
292         !
293         DO jk = 1, jpkm1                 !-- Tracer advective trend
294            DO jj = jplbj+1, jpj-1     
295               DO ji = jplbi+1, jpi-1   ! vector opt.
296                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
297                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
298                  &                                   * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
299               END DO
300           END DO
301         END DO       
302         !                                ! trend diagnostics
303         IF( l_trd )  THEN
304            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
305            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
306         END IF
307         !                                 ! "Poleward" heat and salt transports
308         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
309         !                                 !  heat transport
310         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
311         !
312         !                          !* Vertical advective fluxes
313         !
314         !                                !-- first guess of the slopes
315         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
316         zwx(:,:,jpk) = 0._wp
317         DO jk = 2, jpkm1                       ! interior values
318            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
319         END DO
320         !                                !-- Slopes of tracer
321         zslpx(:,:,1) = 0._wp                   ! surface values
322         DO jk = 2, jpkm1                       ! interior value
323            DO jj = jplbj, jpj
324               DO ji = jplbi, jpi
325                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
326                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
327               END DO
328            END DO
329         END DO
330         DO jk = 2, jpkm1                 !-- Slopes limitation
331            DO jj = jplbj, jpj                      ! interior values
332               DO ji = jplbi, jpi
333                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
334                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
335                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
336               END DO
337            END DO
338         END DO
339         DO jk = 1, jpk-2                 !-- vertical advective flux
340            DO jj = jplbj+1, jpj-1     
341               DO ji = jplbi+1, jpi-1   ! vector opt.
342                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
343                  zalpha = 0.5 + z0w
344                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1)
345                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
346                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
347                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
348               END DO
349            END DO
350         END DO
351         IF( ln_linssh ) THEN                   ! top values, linear free surface only
352            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
353            DO jj = jplbj, jpj     
354               DO ji = jplbi, jpi
355                     zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)
356                  END DO
357               END DO   
358            ELSE                                      ! no cavities: only at the ocean surface
359               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
360            ENDIF
361         ENDIF
362         !
363         DO jk = 1, jpkm1                 !-- vertical advective trend
364            DO jj = jplbj+1, jpj-1     
365               DO ji = jplbi+1, jpi-1   ! vector opt.
366                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
367               END DO
368            END DO
369         END DO
370         !                                ! send trends for diagnostic
371         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
372         !
373      END DO                     ! end of tracer loop
374      !
375#     undef pun
376#     undef pvn
377#     undef pwn
378#     undef ptb
379#     undef pta
380#     undef r1_e1e2t
381#     undef r1_e1e2u
382#     undef r1_e1e2v
383#     undef tmask
384#     undef wmask
385#     undef umask
386#     undef vmask
387#     undef e3u_n
388#     undef e3v_n
389#     undef e3t_n
390#     undef e3w_n
391#     undef mikt
392#     undef rnfmsk
393#     undef upsmsk
394
395   CALL halo_mng_copy(pta_exh2, pta)
396
397    last_khls = jphls - ((SIZE(pta_exh2, 1) - SIZE(pta, 1))/2)
398
399    CALL halo_mng_set(last_khls)
400
401   CALL lbc_lnk( 'traadv_mus', pta, 'T', 1. )
402
403   IF( kt==nitend ) THEN
404      if (allocated(pun_exh2)) DEALLOCATE(pun_exh2)
405      if (allocated(pvn_exh2)) DEALLOCATE(pvn_exh2)
406      if (allocated(pwn_exh2)) DEALLOCATE(pwn_exh2)
407         if (allocated(ptb_exh2)) DEALLOCATE(ptb_exh2)
408      if (allocated(pta_exh2)) DEALLOCATE(pta_exh2)
409      if (allocated(r1_e1e2t_exh2)) DEALLOCATE(r1_e1e2t_exh2)
410      if (allocated(r1_e1e2u_exh2)) DEALLOCATE(r1_e1e2u_exh2)
411      if (allocated(r1_e1e2v_exh2)) DEALLOCATE(r1_e1e2v_exh2)
412      if (allocated(tmask_exh2)) DEALLOCATE(tmask_exh2)
413      if (allocated(wmask_exh2)) DEALLOCATE(wmask_exh2)
414      if (allocated(umask_exh2)) DEALLOCATE(umask_exh2)
415      if (allocated(vmask_exh2)) DEALLOCATE(vmask_exh2)
416      if (allocated(e3u_n_exh2)) DEALLOCATE(e3u_n_exh2)
417      if (allocated(e3v_n_exh2)) DEALLOCATE(e3v_n_exh2)
418      if (allocated(e3t_n_exh2)) DEALLOCATE(e3t_n_exh2)
419      if (allocated(e3w_n_exh2)) DEALLOCATE(e3w_n_exh2)
420         IF (ln_isfcav.and.allocated(mikt_exh2)) DEALLOCATE(mikt_exh2)
421         IF( ld_msc_ups.and.allocated(rnfmsk_exh2)) DEALLOCATE(rnfmsk_exh2)
422      IF( ld_msc_ups.and.allocated(upsmsk_exh2)) DEALLOCATE(upsmsk_exh2)
423
424    ENDIF
425
426   DEALLOCATE(zwx,zwy)
427    DEALLOCATE(zslpx,zslpy)
428   END SUBROUTINE tra_adv_mus
429
430   !!======================================================================
431END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.