source: trunk/NEMO/TOP_SRC/TRP/trcldf_bilap.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: 10.2 KB
Line 
1MODULE trcldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  trcldf_bilap  ***
4   !! TOP :  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if defined key_top
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!----------------------------------------------------------------------
10   !!   trc_ldf_bilap : update the tracer trend with the horizontal diffusion
11   !!                   using a iso-level biharmonic operator
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce_trc         ! ocean dynamics and active tracers variables
15   USE trc             ! ocean passive tracers variables
16   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
17   USE prtctl_trc      ! Print control for debbuging
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC trc_ldf_bilap   ! 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_bilap.F90,v 1.12 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_bilap( kt )
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE trc_ldf_bilap  ***
38      !!
39      !! ** Purpose :   Compute the before horizontal tracer tra diffusive
40      !!      trend and add it to the general trend of tracer equation.
41      !!
42      !! ** Method  :   4th order diffusive operator along model level surfaces
43      !!      evaluated using before fields (forward time scheme). The hor.
44      !!      diffusive trends of passive tracer is given by:
45      !!       * s-coordinate, the vertical scale
46      !!      factors e3. are inside the derivatives:
47      !!      Laplacian of trb:
48      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(trb) ]
49      !!                                  + dj-1[ e1v*e3v/e2v dj(trb) ]  }
50      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
51      !!        zlt   = ahtt * zlt
52      !!        call to lbc_lnk
53      !!      Bilaplacian (laplacian of zlt):
54      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
55      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
56      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
57      !!      Laplacian of trb:
58      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(trb) ]
59      !!                              + dj-1[ e1v/e2v dj(trb) ] }
60      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
61      !!        zlt   = ahtt * zlt
62      !!        call to lbc_lnk
63      !!      Bilaplacian (laplacian of zlt):
64      !!         difft = 1/(e1t*e2t) {  di-1[ e2u/e1u di(zlt) ]
65      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
66      !!
67      !!      Add this trend to the general trend tra :
68      !!         tra = tra + difft
69      !!
70      !! ** Action : - Update tra arrays with the before iso-level
71      !!               biharmonic mixing trend.
72      !!             - Save the trends in trtrd ('key_trc_diatrd')
73      !!
74      !! History :
75      !!        !  91-11  (G. Madec)  Original code
76      !!        !  93-03  (M. Guyon)  symetrical conditions
77      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
78      !!        !  96-01  (G. Madec)  statement function for e3
79      !!        !  96-01  (M. Imbard)  mpp exchange
80      !!        !  97-07  (G. Madec)  optimization, and ahtt
81      !!        !  00-05  (MA Foujols) add lbc for tracer trends
82      !!        !  00-10  (MA Foujols E. Kestenare) use passive tracer coefficient
83      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
84      !!   9.0  !  04-03  (C. Ethe )  F90: Free form and module
85      !!----------------------------------------------------------------------
86      !! * Arguments
87      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
88
89      !! * Local declarations
90      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
91      INTEGER ::   iku, ikv                   ! temporary integers
92
93      REAL(wp) ::   ztra     ! temporary scalars
94
95      REAL(wp), DIMENSION(jpi,jpj) ::   & 
96         zeeu, zeev, zbtr, zlt                 ! workspace
97      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
98         ztu, ztv                              ! workspace
99      CHARACTER (len=22) :: charout
100      !!----------------------------------------------------------------------
101
102      IF( kt == nittrc000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'trc_ldf_bilap : iso-level biharmonic operator'
105         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
106      ENDIF
107      !
108
109      DO jn = 1, jptra
110                                                          ! ===============
111         DO jk = 1, jpkm1                                 ! Horizontal slab
112            !                                             ! ===============
113
114            ! 0. Initialization of metric arrays (for z- or s-coordinates)
115            ! ----------------------------------
116
117            DO jj = 1, jpjm1
118               DO ji = 1, fs_jpim1   ! vector opt.
119#if ! defined key_zco
120                  ! s-coordinates, vertical scale factor are used
121                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
122                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
123                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
124#else
125                  ! z-coordinates, no vertical scale factors
126                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
127                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
128                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
129#endif
130               END DO
131            END DO
132
133
134            ! 1. Laplacian
135            ! ------------
136
137            ! First derivative (gradient)
138            DO jj = 1, jpjm1
139               DO ji = 1, fs_jpim1   ! vector opt.
140                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( trb(ji+1,jj  ,jk,jn) - trb(ji,jj,jk,jn) )
141                  ztv(ji,jj,jk) = zeev(ji,jj) * ( trb(ji  ,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
142               END DO
143            END DO
144
145            IF( ln_zps ) THEN
146               DO jj = 1, jpj-1
147                  DO ji = 1, jpi-1
148                     ! last level
149                     iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
150                     ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
151                     IF( iku == jk ) THEN
152                        ztu(ji,jj,jk) = zeeu(ji,jj) * gtru(ji,jj,jn)
153                     ENDIF
154                     IF( ikv == jk ) THEN
155                        ztv(ji,jj,jk) = zeev(ji,jj) * gtrv(ji,jj,jn)
156                     ENDIF
157                  END DO
158               END DO
159            ENDIF
160
161            ! Second derivative (divergence)
162            DO jj = 2, jpjm1
163               DO ji = fs_2, fs_jpim1   ! vector opt.
164                  zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
165               END DO
166            END DO
167
168            ! Multiply by the eddy diffusivity coefficient
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.
171                  zlt(ji,jj) = fsahtrt(ji,jj,jk) * zlt(ji,jj)
172               END DO
173            END DO
174
175            ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
176            CALL lbc_lnk( zlt, 'T', 1. ) 
177
178            ! 2. Bilaplacian
179            ! --------------
180
181            ! third derivative (gradient)
182            DO jj = 1, jpjm1
183               DO ji = 1, fs_jpim1   ! vector opt.
184                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
185                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
186               END DO
187            END DO
188
189            ! fourth derivative (divergence) and add to the general tracer trend
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1   ! vector opt.
192                  ! horizontal diffusive trends
193                  ztra = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
194                  ! add it to the general tracer trends
195                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
196#if defined key_trc_diatrd
197                  ! save the horizontal diffusive trends
198                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),4) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr(ji,jj)
199                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),5) = (  ztv(ji,jj,jk) - ztv(ji-1,jj,jk) ) * zbtr(ji,jj)
200#endif
201               END DO
202            END DO
203            !                                             ! ===============
204         END DO                                           ! Horizontal slab
205         !                                                ! ===============
206#if defined key_trc_diatrd
207         ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
208         IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),5), 'T', 1. ) 
209#endif
210      END DO
211
212     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
213         WRITE(charout, FMT="('ldf - bilap')")
214         CALL prt_ctl_trc_info(charout)
215         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
216      ENDIF
217
218   END SUBROUTINE trc_ldf_bilap
219
220#else
221   !!----------------------------------------------------------------------
222   !!   Default option                                         Empty module
223   !!----------------------------------------------------------------------
224CONTAINS
225   SUBROUTINE trc_ldf_bilap( kt ) 
226      INTEGER, INTENT(in) :: kt
227      WRITE(*,*) 'trc_ldf_bilap: You should not have seen this print! error?', kt
228   END SUBROUTINE trc_ldf_bilap
229#endif
230   !!==============================================================================
231END MODULE trcldf_bilap
Note: See TracBrowser for help on using the repository browser.