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.
traldf_iso.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90 @ 4431

Last change on this file since 4431 was 4427, checked in by trackstand2, 10 years ago

First files changed on last FINISS work package. Stephen's work although
commited by Andy P.
  • Property svn:keywords set to Id
File size: 25.3 KB
Line 
1MODULE traldf_iso
2   !!======================================================================
3   !!                   ***  MODULE  traldf_iso  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :  OPA  !  1994-08  (G. Madec, M. Imbard)
7   !!            8.0  !  1997-05  (G. Madec)  split into traldf and trazdf
8   !!            NEMO !  2002-08  (G. Madec)  Free form, F90
9   !!            1.0  !  2005-11  (G. Madec)  merge traldf and trazdf :-)
10   !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
11   !!----------------------------------------------------------------------
12#if   defined key_ldfslp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_ldfslp'               slope of the lateral diffusive direction
15   !!----------------------------------------------------------------------
16   !!   tra_ldf_iso  : update the tracer trend with the horizontal
17   !!                  component of a iso-neutral laplacian operator
18   !!                  and with the vertical part of
19   !!                  the isopycnal or geopotential s-coord. operator
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and active tracers
22   USE dom_oce         ! ocean space and time domain
23   USE trc_oce         ! share passive tracers/Ocean variables
24   USE zdf_oce         ! ocean vertical physics
25   USE ldftra_oce      ! ocean active tracers: lateral physics
26   USE ldfslp          ! iso-neutral slopes
27   USE diaptr          ! poleward transport diagnostics
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O library
30#if defined key_diaar5
31   USE phycst          ! physical constants
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33#endif
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_ldf_iso   ! routine called by step.F90
39
40   !! * Control permutation of array indices
41#  include "oce_ftrans.h90"
42#  include "dom_oce_ftrans.h90"
43#  include "trc_oce_ftrans.h90"
44#  include "zdf_oce_ftrans.h90"
45#  include "ldftra_oce_ftrans.h90"
46#  include "ldfslp_ftrans.h90"
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "ldftra_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv,              &
60      &                                ptb, pta, kjpt, pahtb0 )
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE tra_ldf_iso  ***
63      !!
64      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
65      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
66      !!      add it to the general trend of tracer equation.
67      !!
68      !! ** Method  :   The horizontal component of the lateral diffusive trends
69      !!      is provided by a 2nd order operator rotated along neural or geopo-
70      !!      tential surfaces to which an eddy induced advection can be added
71      !!      It is computed using before fields (forward in time) and isopyc-
72      !!      nal or geopotential slopes computed in routine ldfslp.
73      !!
74      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
75      !!      ========    with partial cell update if ln_zps=T.
76      !!
77      !!      2nd part :  horizontal fluxes of the lateral mixing operator
78      !!      ========   
79      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
80      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
81      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
82      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
83      !!      take the horizontal divergence of the fluxes:
84      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
85      !!      Add this trend to the general trend (ta,sa):
86      !!         ta = ta + difft
87      !!
88      !!      3rd part: vertical trends of the lateral mixing operator
89      !!      ========  (excluding the vertical flux proportional to dk[t] )
90      !!      vertical fluxes associated with the rotated lateral mixing:
91      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
92      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
93      !!      take the horizontal divergence of the fluxes:
94      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
95      !!      Add this trend to the general trend (ta,sa):
96      !!         pta = pta + difft
97      !!
98      !! ** Action :   Update pta arrays with the before rotated diffusion
99      !!----------------------------------------------------------------------
100!      USE arpdebugging, ONLY: dump_array
101      USE timing,   ONLY: timing_start, timing_stop
102      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
103      USE oce     , ONLY:   zftu => ua       , zftv  => va         ! (ua,va) used as workspace
104      !! DCSE_NEMO: need additional directives for renamed module variables
105!FTRANS zftu zftv :I :I :z
106#if defined key_z_first
107!      USE wrk_nemo, ONLY:   wdkt => wrk_3d_9 , wdk1t => wrk_3d_10  ! 3D workspace
108!FTRANS wdkt wdk1t :I :I :z
109#else
110      USE wrk_nemo, ONLY:   zdkt => wrk_2d_1 , zdk1t => wrk_2d_2
111#endif
112      USE wrk_nemo, ONLY:   z2d  => wrk_2d_3                       ! 2D workspace
113      USE wrk_nemo, ONLY:   zdit => wrk_3d_6 , zdjt  => wrk_3d_7 , ztfw => wrk_3d_8   ! 3D workspace
114!FTRANS zdit zdjt ztfw :I :I :z
115
116      !
117      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
118      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
119      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
120      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
121
122      !! DCSE_NEMO: This style defeats ftrans
123!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
124!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
125!FTRANS ptb pta :I :I :z :
126      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)        ! before and now tracer fields
127      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)        ! tracer trend
128
129      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
130      !
131      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
132      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
133      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
134      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      -
135#if defined key_z_first
136      REAL(wp) ::  wdkt , wdki1t , wdkim1t , wdkj1t , wdkjm1t
137      REAL(wp) ::  wdk1t, wdk1i1t, wdk1im1t, wdk1j1t, wdk1jm1t
138#endif
139
140#if defined key_diaar5
141      REAL(wp)                         ::   zztmp               ! local scalar
142#endif
143      !!----------------------------------------------------------------------
144
145      CALL timing_start('tra_ldf_iso')
146
147#if defined key_z_first
148!! SMP Should not need to reserve or release 9 and 10 any more.
149      IF( wrk_in_use(3, 6,7,8,9,10) .OR. wrk_in_use(2, 3) ) THEN
150#else
151      IF( wrk_in_use(3, 6,7,8) .OR. wrk_in_use(2, 1,2,3) ) THEN
152#endif
153          CALL ctl_stop('tra_ldf_iso : requested workspace array unavailable')   ;   RETURN
154      ENDIF
155
156      IF( kt == nit000 )  THEN
157         IF(lwp) WRITE(numout,*)
158         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
159         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
160      ENDIF
161      !
162!      CALL dump_array(kt, 'ptb', ptb(:,:,1,1), withHalos=.TRUE.)
163      !                                                          ! ===========
164!DIR$ SHORTLOOP
165
166#if defined key_z_first
167      zdit(:,:,:) = 0.0_wp
168      zdjt(:,:,:) = 0.0_wp
169#endif
170
171      DO jn = 1, kjpt                                            ! tracer loop
172         !                                                       ! ===========
173         !                                               
174         !!----------------------------------------------------------------------
175         !!   I - masked horizontal derivative
176         !!----------------------------------------------------------------------
177         !CALL timing_start('traldf_iso_I')
178
179         ! Horizontal tracer gradient
180#if defined key_z_first
181         DO jj = 1, jpjm1
182            DO ji = 1, jpim1
183               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
184                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
185                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
186               END DO
187            END DO
188         END DO
189#else
190         !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
191         zdit (1,:,1:jpkf) = 0.e0     ;     zdit (jpi,:,1:jpkf) = 0.e0
192         zdjt (1,:,1:jpkf) = 0.e0     ;     zdjt (jpi,:,1:jpkf) = 0.e0
193         !!end
194         DO jk = 1, jpkfm1 ! jpkm1
195            DO jj = 1, jpjm1
196               DO ji = 1, fs_jpim1   ! vector opt.
197                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
198                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
199               END DO
200            END DO
201         END DO
202#endif
203         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level
204            DO jj = 1, jpjm1
205               DO ji = 1, fs_jpim1   ! vector opt.
206                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
207                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)     
208               END DO
209            END DO
210         ENDIF
211         !
212         !CALL timing_stop('traldf_iso_I','section')
213
214         !!----------------------------------------------------------------------
215         !!   II - horizontal trend  (full)
216         !!----------------------------------------------------------------------
217         !CALL timing_start('traldf_iso_II')
218#if defined key_z_first
219            ! 1. Vertical tracer gradient at level jk and jk+1
220            ! ------------------------------------------------
221            ! surface boundary condition: wdkt(jk=1)=wdkt(jk=2)
222
223!!$         DO jj = 1, jpj
224!!$            DO ji = 1, jpi
225!!$               DO jk = 1, jpkm1
226!!$                  wdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
227!!$               END DO
228!!$               wdkt(ji,jj,1) = wdk1t(ji,jj,1)
229!!$               DO jk = 2, jpkm1
230!!$                  wdkt(ji,jj,jk) =  ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
231!!$               END DO
232!!$            END DO
233!!$         END DO
234
235            ! 2. Horizontal fluxes
236            ! --------------------   
237!!$         DO jj = 1 , jpjm1
238!!$            DO ji = 1, jpim1
239!!$               DO jk = 1, jpkm1
240!!$                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
241!!$                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
242!!$                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
243!!$                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
244!!$                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
245!!$                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
246!!$                  zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
247!!$                  zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
248!!$                  zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk)   &
249!!$                     &              + zcof1 * (  wdkt (ji+1,jj,jk) + wdk1t(ji,jj,jk)      &
250!!$                     &                         + wdk1t(ji+1,jj,jk) + wdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk)
251!!$                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
252!!$                     &              + zcof2 * (  wdkt (ji,jj+1,jk) + wdk1t(ji,jj,jk)      &
253!!$                     &                         + wdk1t(ji,jj+1,jk) + wdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                 
254!!$               END DO
255!!$            END DO
256!!$         END DO
257
258         DO jj = 2 , jpjm1
259            DO ji = 2, jpim1
260               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
261
262                  ! 1. Vertical tracer gradient at level jk and jk+1
263                  ! ------------------------------------------------
264                  ! surface boundary condition: wdkt(jk=1)=wdkt(jk=2)
265
266                  wdk1t = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1)
267                  wdk1i1t = ( ptb(ji+1,jj,jk,jn) - ptb(ji+1,jj,jk+1,jn) ) * tmask(ji+1,jj,jk+1)
268                  wdk1im1t = ( ptb(ji-1,jj,jk,jn) - ptb(ji-1,jj,jk+1,jn) ) * tmask(ji-1,jj,jk+1)
269                  wdk1j1t = ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj+1,jk+1,jn) ) * tmask(ji,jj+1,jk+1)
270                  wdk1jm1t = ( ptb(ji,jj-1,jk,jn) - ptb(ji,jj-1,jk+1,jn) ) * tmask(ji,jj-1,jk+1)
271
272                  IF(jk > 1)THEN
273                     wdkt =  ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
274                     wdki1t = ( ptb(ji+1,jj,jk-1,jn) - ptb(ji+1,jj,jk,jn) ) * tmask(ji+1,jj,jk)
275                     wdkim1t = ( ptb(ji-1,jj,jk-1,jn) - ptb(ji-1,jj,jk,jn) ) * tmask(ji-1,jj,jk)
276                     wdkj1t = ( ptb(ji,jj+1,jk-1,jn) - ptb(ji,jj+1,jk,jn) ) * tmask(ji,jj+1,jk)
277                     wdkjm1t = ( ptb(ji,jj-1,jk-1,jn) - ptb(ji,jj-1,jk,jn) ) * tmask(ji,jj-1,jk)
278                  ELSE
279                     wdkt   = wdk1t
280                     wdki1t = wdk1i1t
281                     wdkim1t= wdk1im1t
282                     wdkj1t = wdk1j1t
283                     wdkjm1t= wdk1jm1t
284                  END IF
285
286                  ! II.4 Second derivative (divergence) and add to the general trend
287                  ! ----------------------------------------------------------------
288                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
289
290                  ztra = zbtr * (                                                   &
291
292!                                zftu(ji,jj,jk) -
293                            ( ((fsahtu(ji,jj,jk) + pahtb0) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)) * zdit(ji,jj,jk)                                &
294                              - ( fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) /                     &
295                                MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)                          &
296                                    + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1.) ) *                 &
297                              (wdki1t + wdk1t + wdk1i1t + wdkt) ) * umask(ji,jj,jk) - &
298
299!                                zftu(ji-1,jj,jk) +
300                            ( ((fsahtu(ji-1,jj,jk) + pahtb0) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) / e1u(ji-1,jj)) * zdit(ji-1,jj,jk)                                &
301                              - ( fsahtu(ji-1,jj,jk) * e2u(ji-1,jj) * uslp(ji-1,jj,jk) /                 &
302                                MAX(  tmask(ji,jj,jk  ) + tmask(ji-1,jj,jk+1)                            &
303                                    + tmask(ji,jj,jk+1) + tmask(ji-1,jj,jk  ), 1.) ) *                   &
304                              (wdkt + wdk1im1t + wdk1t + wdkim1t) ) * umask(ji-1,jj,jk) + &
305
306
307!                                zftv(ji,jj,jk) -
308                            (  ((fsahtv(ji,jj,jk) + pahtb0) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)) * zdjt(ji,jj,jk)   &
309                     &              - ( fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) /                 &
310                                MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)                            &
311                     &              + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )) *                   &
312                                (wdkj1t + wdk1t + wdk1j1t + wdkt)  ) * vmask(ji,jj,jk) - &
313!                                zftv(ji,jj-1,jk) &
314                            (  ((fsahtv(ji,jj-1,jk) + pahtb0) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) / e2v(ji,jj-1)) * zdjt(ji,jj-1,jk)   &
315                     &              - ( fsahtv(ji,jj-1,jk) * e1v(ji,jj-1) * vslp(ji,jj-1,jk) /           &
316                                MAX(  tmask(ji,jj,jk  ) + tmask(ji,jj-1,jk+1)                            &
317                     &              + tmask(ji,jj,jk+1) + tmask(ji,jj-1,jk  ), 1. )) *                   &
318                                (wdkt + wdk1jm1t + wdk1t + wdkjm1t)  ) * vmask(ji,jj-1,jk) &
319
320                                )
321                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
322               END DO
323            END DO
324         END DO
325#else
326!CDIR PARALLEL DO PRIVATE( zdk1t )
327         !                                                ! ===============
328         DO jk = 1, jpkfm1 ! jpkm1                        ! Horizontal slab
329            !                                             ! ===============
330            ! 1. Vertical tracer gradient at level jk and jk+1
331            ! ------------------------------------------------
332            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
333            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
334            !
335            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)
336            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
337            ENDIF
338
339            ! 2. Horizontal fluxes
340            ! --------------------   
341            DO jj = 1 , jpjm1
342               DO ji = 1, fs_jpim1   ! vector opt.
343                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
344                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
345                  !
346                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
347                     &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
348                  !
349                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
350                     &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
351                  !
352                  zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
353                  zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
354                  !
355                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
356                     &              + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
357                     &                         + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
358                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
359                     &              + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
360                     &                         + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                 
361               END DO
362            END DO
363
364            ! II.4 Second derivative (divergence) and add to the general trend
365            ! ----------------------------------------------------------------
366            DO jj = 2 , jpjm1
367               DO ji = fs_2, fs_jpim1   ! vector opt.
368                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
369                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
370                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
371               END DO
372            END DO
373            !                                          ! ===============
374         END DO                                        !   End of slab 
375         !                                             ! ===============
376#endif
377         !
378         ! "Poleward" diffusive heat or salt transports (T-S case only)
379         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
380            IF( jn == jp_tem)THEN
381               htr_ldf = ptr_vj( zftv )
382            END IF
383            IF( jn == jp_sal)THEN
384               str_ldf = ptr_vj( zftv )
385            END IF
386         ENDIF
387 
388#if defined key_diaar5
389         IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
390            z2d(:,:) = 0._wp 
391            zztmp = rau0 * rcp 
392#if defined key_z_first
393            DO jj = 2, jpjm1
394               DO ji = 2, jpim1
395                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
396#else
397            DO jk = 1, jpkfm1 ! jpkm1
398               DO jj = 2, jpjm1
399                  DO ji = fs_2, fs_jpim1   ! vector opt.
400#endif
401                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
402                  END DO
403               END DO
404            END DO
405            z2d(:,:) = zztmp * z2d(:,:)
406            CALL lbc_lnk( z2d, 'U', -1. )
407            CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
408            z2d(:,:) = 0._wp 
409#if defined key_z_first
410            DO jj = 2, jpjm1
411               DO ji = 2, jpim1
412                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
413#else
414            DO jk = 1, jpkfm1 ! jpkm1
415               DO jj = 2, jpjm1
416                  DO ji = fs_2, fs_jpim1   ! vector opt.
417#endif
418                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
419                  END DO
420               END DO
421            END DO
422            z2d(:,:) = zztmp * z2d(:,:)
423            CALL lbc_lnk( z2d, 'V', -1. )
424            CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
425         END IF
426#endif
427         !CALL timing_stop('traldf_iso_II','section')
428
429         !!--------------------------------------------------------------------
430         !!   III - vertical trend of T & S (extra diagonal terms only)
431         !!--------------------------------------------------------------------
432         !CALL timing_start('traldf_iso_III')
433         
434         ! Local constant initialization
435         ! -----------------------------
436         ! Vertical fluxes
437         ! ---------------
438         
439         ! Surface and bottom vertical fluxes set to zero
440#if defined key_z_first
441         ztfw(:,:,:) = 0.0_wp
442#else
443         ztfw(1,:,1:jpkf) = 0.e0     ;     ztfw(jpi,:,1:jpkf) = 0.e0
444         ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpkf) = 0.e0
445#endif
446
447         ! interior (2=<jk=<jpk-1)
448#if defined key_z_first
449         DO jj = 2, jpjm1
450            DO ji = 2, jpim1
451               DO jk = 2, mbkmax(ji,jj)-1
452#else
453         DO jk = 2, jpkfm1
454            DO jj = 2, jpjm1
455               DO ji = fs_2, fs_jpim1   ! vector opt.
456#endif
457                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
458                  !
459                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
460                     &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
461                  zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
462                     &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
463                  !
464                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)
465                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
466                  !
467                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
468                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
469                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
470                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
471               END DO
472            END DO
473         END DO
474         
475         
476         ! I.5 Divergence of vertical fluxes added to the general tracer trend
477         ! -------------------------------------------------------------------
478#if defined key_z_first
479         DO jj = 2, jpjm1
480            DO ji = 2, jpim1
481               DO jk = 1, mbkmax(ji,jj)-1
482#else
483         DO jk = 1, jpkfm1
484            DO jj = 2, jpjm1
485               DO ji = fs_2, fs_jpim1   ! vector opt.
486#endif
487                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
488                  ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
489                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
490               END DO
491            END DO
492         END DO
493         !
494
495         !CALL timing_stop('traldf_iso_III','section')
496
497      END DO
498      !
499#if defined key_z_first
500      IF( wrk_not_released(3, 6,7,8,9,10) .OR.   &
501          wrk_not_released(2, 3) )       CALL ctl_stop('tra_ldf_iso: failed to release workspace arrays')
502#else
503      IF( wrk_not_released(3, 6,7,8) .OR.   &
504          wrk_not_released(2, 1,2,3) )   CALL ctl_stop('tra_ldf_iso: failed to release workspace arrays')
505#endif
506      !
507      CALL timing_stop('tra_ldf_iso','section')
508      !
509   END SUBROUTINE tra_ldf_iso
510
511#else
512   !!----------------------------------------------------------------------
513   !!   default option :   Dummy code   NO rotation of the diffusive tensor
514   !!----------------------------------------------------------------------
515CONTAINS
516   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 )      ! Empty routine
517      CHARACTER(len=3) ::   cdtype
518      REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels
519      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
520      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, cdtype, pgu(1,1,1), pgv(1,1,1),   &
521         &                                                             ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0
522   END SUBROUTINE tra_ldf_iso
523#endif
524
525   !!==============================================================================
526END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.