source: trunk/NEMO/TOP_SRC/TRP/trcldf_iso_zps.F90 @ 941

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

  • Property svn:executable set to *
File size: 12.3 KB
Line 
1MODULE trcldf_iso_zps
2   !!==============================================================================
3   !!                   ***  MODULE  trcldf_iso_zps  ***
4   !! Ocean passive tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if defined key_top &&  defined key_ldfslp 
7   !!----------------------------------------------------------------------
8   !!   'key_top'          and                                   TOP models
9   !!   'key_ldfslp'               slope of the lateral diffusive direction
10   !!----------------------------------------------------------------------
11   !!   trc_ldf_iso_zps : update the tracer trend with the horizontal
12   !!                     component of a iso-neutral laplacian operator
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce_trc             ! ocean dynamics and active tracers variables
16   USE trc                 ! ocean passive tracers variables
17   USE prtctl_trc          ! Print control for debbuging
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Accessibility
23   PUBLIC trc_ldf_iso_zps  ! routine called by step.F90
24
25   !! * Substitutions
26#  include "top_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   TOP 1.0 , LOCEAN-IPSL (2005)
29   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcldf_iso_zps.F90,v 1.11 2007/10/12 09:26:30 opalod Exp $
30   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
31   !!----------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE trc_ldf_iso_zps( kt )
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE trc_ldf_iso_zps  ***
38      !!
39      !! ** Purpose :   Compute the before horizontal tracer  diffusive
40      !!      trend and add it to the general trend of tracer equation.
41      !!
42      !! ** Method  :   The horizontal component of the lateral diffusive trends
43      !!      is provided by a 2nd order operator rotated along neural or geopo-
44      !!      tential surfaces to which an eddy induced advection can be added
45      !!      It is computed using before fields (forward in time) and isopyc-
46      !!      nal or geopotential slopes computed in routine ldfslp.
47      !!
48      !!      horizontal fluxes associated with the rotated lateral mixing:
49      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
50      !!               - aht       e2u*uslp    dk[ mi(mk(trb)) ]
51      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
52      !!               - aht       e2u*vslp    dk[ mj(mk(trb)) ]
53      !!      add horizontal Eddy Induced advective fluxes (lk_traldf_eiv=T):
54      !!         zftu = zftu - dk-1[ aht e2u mi(wslpi) ] mi( trb )
55      !!         zftv = zftv - dk-1[ aht e1v mj(wslpj) ] mj( trb )
56      !!      take the horizontal divergence of the fluxes:
57      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
58      !!      Add this trend to the general trend tra :
59      !!         tra = tra + difft
60      !!
61      !!      'key_trdtra' defined: the trend is saved for diagnostics.
62      !!
63      !!      macro-tasked on horizontal slab (jk-loop).
64      !!
65      !! ** Action :
66      !!         Update tra arrays with the before along level biharmonic
67      !!      mixing trend.
68      !!         Save in trtrd arrays the trends if 'key_trc_diatrd' defined
69      !!
70      !! History :
71      !!        !  94-08  (G. Madec, M. Imbard)
72      !!        !  97-05  (G. Madec)  split into traldf and trazdf
73      !!   8.5  !  02-08  (G. Madec)  Free form, F90
74      !!   9.0  !  04-03  (C. Ethe)  adapted for passive tracers
75      !!----------------------------------------------------------------------
76      !! * Modules used
77      USE oce_trc       , zftu => ua,  &  ! use ua as workspace
78         &                zfsu => va      ! use va as workspace
79
80      !! * Arguments
81      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
82
83      !! * Local declarations
84      INTEGER ::   ji, jj, jk,jn          ! dummy loop indices
85      INTEGER ::   iku, ikv               ! temporary integer
86      REAL(wp) ::   &
87         zabe1, zabe2, zcof1, zcof2,   &  ! temporary scalars
88         zmsku, zmskv, zbtr, ztra       
89
90      REAL(wp), DIMENSION(jpi,jpj) ::   &
91         zdkt , zdk1t                     ! temporary workspace
92
93      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
94         zftv, zgtbu, zgtbv               ! temporary workspace
95
96#if defined key_trcldf_eiv
97      REAL(wp), DIMENSION(jpi,jpj) ::   &
98         zftug, zftvg                    ! temporary workspace
99
100      REAL(wp) ::   &
101         zuwk, zvwk,                   &
102         zuwk1, zvwk1,                 &
103         zcg1,zcg2
104
105#if defined key_trc_diatrd
106      REAL(wp) ::   &
107         ztagu, ztagv
108#endif
109
110#endif
111      CHARACTER (len=22) :: charout
112      !!----------------------------------------------------------------------
113
114      IF( kt == nittrc000 ) THEN
115         IF(lwp) WRITE(numout,*)
116         IF(lwp) WRITE(numout,*) 'trc_ldf_iso_zps : iso neutral laplacian diffusion in '
117         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   z-coordinates with partial steps'
118#if defined key_trcldf_eiv && defined key_diaeiv
119         u_trc_eiv(:,:,:) = 0.e0
120         v_trc_eiv(:,:,:) = 0.e0
121#endif
122      ENDIF
123
124      DO jn = 1, jptra
125
126#if defined key_trcldf_eiv && key_trc_diatrd
127         ztagu = 0.e0
128         ztagv = 0.e0
129#endif
130
131         ! Horizontal passive tracer gradient
132         DO jk = 1, jpk
133            DO jj = 1, jpj-1
134               DO ji = 1, fs_jpim1   ! vector opt.
135                  zgtbu(ji,jj,jk) = tmask(ji,jj,jk) * ( trb(ji+1,jj  ,jk,jn) - trb(ji,jj,jk,jn) )
136                  zgtbv(ji,jj,jk) = tmask(ji,jj,jk) * ( trb(ji  ,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
137               END DO
138            END DO
139         END DO
140         ! partial steps correction at the last level
141         DO jj = 1, jpj-1
142            DO ji = 1, jpi-1
143               ! last level
144               iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
145               ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
146               zgtbu(ji,jj,iku) = gtru(ji,jj,jn)               
147               zgtbv(ji,jj,ikv) = gtrv(ji,jj,jn)               
148            END DO
149         END DO
150
151         !                                                ! ===============
152         DO jk = 1, jpkm1                                 ! Horizontal slab
153            !                                             ! ===============
154            ! 1. Vertical tracer gradient at level jk and jk+1
155            ! ------------------------------------------------
156            ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
157
158            zdk1t(:,:) = ( trb(:,:,jk,jn) - trb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
159
160            IF( jk == 1 ) THEN
161               zdkt(:,:) = zdk1t(:,:)
162            ELSE
163               zdkt(:,:) = ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) ) * tmask(:,:,jk)
164            ENDIF
165
166
167            ! 2. Horizontal fluxes
168            ! --------------------
169
170            DO jj = 1 , jpjm1
171               DO ji = 1, fs_jpim1   ! vector opt.
172                  zabe1 = ( fsahtru(ji,jj,jk) + ahtrb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
173                  zabe2 = ( fsahtrv(ji,jj,jk) + ahtrb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
174
175                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
176                     + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
177
178                  zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
179                     + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
180
181                  zcof1 = -fsahtru(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
182                  zcof2 = -fsahtrv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
183
184                  zftu(ji,jj,jk) = umask(ji,jj,jk) * (  zabe1 * zgtbu(ji,jj,jk)   &
185                     &                                + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
186                     &                                           + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  )
187                  zftv(ji,jj,jk) = vmask(ji,jj,jk) * (  zabe2 * zgtbv(ji,jj,jk)   &
188                     &                                + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
189                     &                                           + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  )
190               END DO
191            END DO
192
193# if defined key_trcldf_eiv
194            ! ---------------------------------------!
195            ! Eddy induced vertical advective fluxes !
196            ! ---------------------------------------!
197            DO jj = 1, jpjm1
198               DO ji = 1, fs_jpim1   ! vector opt.
199                  zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj  ,jk  ) ) * fsaeitru(ji,jj,jk  ) * umask(ji,jj,jk  )
200                  zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj  ,jk+1) ) * fsaeitru(ji,jj,jk+1) * umask(ji,jj,jk+1)
201                  zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji  ,jj+1,jk  ) ) * fsaeitrv(ji,jj,jk  ) * vmask(ji,jj,jk  )
202                  zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji  ,jj+1,jk+1) ) * fsaeitrv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
203
204                  zcg1= -0.25 * e2u(ji,jj) * umask(ji,jj,jk) * ( zuwk-zuwk1 )
205                  zcg2= -0.25 * e1v(ji,jj) * vmask(ji,jj,jk) * ( zvwk-zvwk1 )
206
207                  zftug(ji,jj) = zcg1 * ( trb(ji+1,jj,jk,jn) + trb(ji,jj,jk,jn) )
208                  zftvg(ji,jj) = zcg2 * ( trb(ji,jj+1,jk,jn) + trb(ji,jj,jk,jn) )
209
210                  zftu(ji,jj,jk) = zftu(ji,jj,jk) + zftug(ji,jj)
211                  zftv(ji,jj,jk) = zftv(ji,jj,jk) + zftvg(ji,jj)
212
213#   if defined key_diaeiv
214                  u_trc_eiv(ji,jj,jk) = -2. * zcg1 / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
215                  v_trc_eiv(ji,jj,jk) = -2. * zcg2 / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
216#   endif
217               END DO
218            END DO
219# endif
220
221            ! II.4 Second derivative (divergence) and add to the general trend
222            ! ----------------------------------------------------------------
223
224            DO jj = 2 , jpjm1
225               DO ji = fs_2, fs_jpim1   ! vector opt.
226                  zbtr= 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
227                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
228                  tra (ji,jj,jk,jn) = tra (ji,jj,jk,jn) + ztra
229#if defined key_trc_diatrd
230                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),4) = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,  jj,jk) )
231                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),5) = zbtr * ( zftv(ji,jj,jk) - zftv(ji  ,jj-1,jk) )
232#endif
233               END DO
234            END DO
235#if defined key_trc_diatrd
236#   if defined key_trcldf_eiv
237            DO jj = 2 , jpjm1
238               DO ji = fs_2, fs_jpim1   ! vector opt.
239                  zbtr= 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
240                  ztagu = ( zftug(ji,jj) - zftug(ji-1,jj  ) ) * zbtr
241                  ztagv = ( zftvg(ji,jj) - zftvg(ji  ,jj-1) ) * zbtr
242                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),4) = trtrd(ji,jj,jk,ikeep(jn),4) - ztagu
243                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),5) = trtrd(ji,jj,jk,ikeep(jn),5) - ztagv
244                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),7) = ztagu
245                  IF (luttrd(jn)) trtrd (ji,jj,jk,ikeep(jn),8) = ztagv
246               END DO
247            END DO
248#   endif
249#endif
250            !                                          ! ===============
251         END DO                                        !   End of slab 
252         !                                             ! ===============
253      END DO
254
255      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
256         WRITE(charout, FMT="('ldf - iso/zps')")
257         CALL prt_ctl_trc_info(charout)
258         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
259      ENDIF
260
261   END SUBROUTINE trc_ldf_iso_zps
262
263#else
264   !!----------------------------------------------------------------------
265   !!   default option :   Dummy code   NO rotation of the diffusive tensor
266   !!----------------------------------------------------------------------
267CONTAINS
268   SUBROUTINE trc_ldf_iso_zps( kt )               ! Empty routine
269      INTEGER, INTENT(in) :: kt
270      WRITE(*,*) 'trc_ldf_iso_zps: You should not have seen this print! error?', kt
271   END SUBROUTINE trc_ldf_iso_zps
272#endif
273
274   !!==============================================================================
275END MODULE trcldf_iso_zps
Note: See TracBrowser for help on using the repository browser.