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_bilap.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traldf_bilap.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 KB
Line 
1MODULE traldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  traldf_bilap  ***
4   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   tra_ldf_bilap : update the tracer trend with the horizontal diffusion
9   !!                   using a iso-level biharmonic operator
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and active tracers
13   USE dom_oce         ! ocean space and time domain
14   USE ldftra_oce     ! ocean tracer   lateral physics
15   USE trdtra_oce     ! ocean active tracer trend
16   USE in_out_manager  ! I/O manager
17   USE ldfslp          ! iso-neutral slopes
18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Routine accessibility
24   PUBLIC tra_ldf_bilap   ! routine called by step.F90
25
26   !! * Substitutions
27#  include "domzgr_substitute.h90"
28#  include "ldftra_substitute.h90"
29#  include "ldfeiv_substitute.h90"
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !!   OPA 9.0 , LODYC-IPSL (2003)
33   !!----------------------------------------------------------------------
34
35CONTAINS
36   
37   SUBROUTINE tra_ldf_bilap( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE tra_ldf_bilap  ***
40      !!
41      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
42      !!      trend and add it to the general trend of tracer equation.
43      !!
44      !! ** Method  :   4th order diffusive operator along model level surfaces
45      !!      evaluated using before fields (forward time scheme). The hor.
46      !!      diffusive trends of temperature (idem for salinity) is given by:
47      !!       * s-coordinate ('key_s_coord' defined), the vertical scale
48      !!      factors e3. are inside the derivatives:
49      !!      Laplacian of tb:
50      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
51      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
52      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
53      !!        zlt   = ahtt * zlt
54      !!        call to lbc_lnk
55      !!      Bilaplacian (laplacian of zlt):
56      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
57      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
58      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
59      !!      Laplacian of tb:
60      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(tb) ]
61      !!                              + dj-1[ e1v/e2v dj(tb) ] }
62      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
63      !!        zlt   = ahtt * zlt
64      !!        call to lbc_lnk
65      !!      Bilaplacian (laplacian of zlt):
66      !!         difft = 1/(e1t*e2t) {  di-1[ e2u/e1u di(zlt) ]
67      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
68      !!
69      !!      Add this trend to the general trend (ta,sa):
70      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
71      !!
72      !! ** Action : - Update (ta,sa) arrays with the before iso-level
73      !!               biharmonic mixing trend.
74      !!             - Save the trends in (ttrd,strd) ('key_diatrends')
75      !!
76      !! History :
77      !!        !  91-11  (G. Madec)  Original code
78      !!        !  93-03  (M. Guyon)  symetrical conditions
79      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
80      !!        !  96-01  (G. Madec)  statement function for e3
81      !!        !  96-01  (M. Imbard)  mpp exchange
82      !!        !  97-07  (G. Madec)  optimization, and ahtt
83      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
84      !!----------------------------------------------------------------------
85      !! * Arguments
86      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
87
88      !! * Local declarations
89      INTEGER ::   ji, jj, jk             ! dummy loop indices
90#if defined key_partial_steps
91      INTEGER ::   iku, ikv               ! temporary integers
92#endif
93      REAL(wp) ::   zta, zsa    ! temporary scalars
94      REAL(wp), DIMENSION(jpi,jpj) ::   & 
95         zeeu, zeev, zbtr,             &  ! workspace
96         zlt, ztu, ztv,                &
97         zls, zsu, zsv
98      !!----------------------------------------------------------------------
99
100      IF( kt == nit000 ) THEN
101         IF(lwp) WRITE(numout,*)
102         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
103         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
104      ENDIF
105      !                                                ! ===============
106      DO jk = 1, jpkm1                                 ! Horizontal slab
107         !                                             ! ===============
108
109         ! 0. Initialization of metric arrays (for z- or s-coordinates)
110         ! ----------------------------------
111
112         DO jj = 1, jpjm1
113            DO ji = 1, fs_jpim1   ! vector opt.
114#if defined key_s_coord || defined key_partial_steps
115               ! s-coordinates, vertical scale factor are used
116               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
117               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
118               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
119#else
120               ! z-coordinates, no vertical scale factors
121               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
122               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
123               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
124#endif
125            END DO
126         END DO
127
128
129         ! 1. Laplacian
130         ! ------------
131
132         ! First derivative (gradient)
133         DO jj = 1, jpjm1
134            DO ji = 1, fs_jpim1   ! vector opt.
135               ztu(ji,jj) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
136               zsu(ji,jj) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
137               ztv(ji,jj) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
138               zsv(ji,jj) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
139            END DO
140         END DO
141#if defined key_partial_steps
142         DO jj = 1, jpj-1
143            DO ji = 1, jpi-1
144               ! last level
145               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
146               ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
147               IF( iku == jk ) THEN
148                  ztu(ji,jj) = zeeu(ji,jj) * gtu(ji,jj)
149                  zsu(ji,jj) = zeeu(ji,jj) * gsu(ji,jj)
150               ENDIF
151               IF( ikv == jk ) THEN
152                  ztv(ji,jj) = zeev(ji,jj) * gtv(ji,jj)
153                  zsv(ji,jj) = zeev(ji,jj) * gsv(ji,jj)
154               ENDIF
155            END DO
156         END DO
157#endif
158
159         ! Second derivative (divergence)
160         DO jj = 2, jpjm1
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj) - ztu(ji-1,jj) + ztv(ji,jj) - ztv(ji,jj-1)  )
163               zls(ji,jj) = zbtr(ji,jj) * (  zsu(ji,jj) - zsu(ji-1,jj) + zsv(ji,jj) - zsv(ji,jj-1)  )
164            END DO
165         END DO
166
167         ! Multiply by the eddy diffusivity coefficient
168         DO jj = 2, jpjm1
169            DO ji = fs_2, fs_jpim1   ! vector opt.
170               zlt(ji,jj) = fsahtt(ji,jj,jk) * zlt(ji,jj)
171               zls(ji,jj) = fsahtt(ji,jj,jk) * zls(ji,jj)
172            END DO
173         END DO
174
175         ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
176         CALL lbc_lnk( zlt, 'T', 1. )   ;    CALL lbc_lnk( zls, '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) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
185               zsu(ji,jj) = zeeu(ji,jj) * ( zls(ji+1,jj  ) - zls(ji,jj) )
186               ztv(ji,jj) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
187               zsv(ji,jj) = zeev(ji,jj) * ( zls(ji  ,jj+1) - zls(ji,jj) )
188            END DO
189         END DO
190
191         ! fourth derivative (divergence) and add to the general tracer trend
192         DO jj = 2, jpjm1
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               ! horizontal diffusive trends
195               zta = zbtr(ji,jj) * (  ztu(ji,jj) - ztu(ji-1,jj) + ztv(ji,jj) - ztv(ji,jj-1)  )
196               zsa = zbtr(ji,jj) * (  zsu(ji,jj) - zsu(ji-1,jj) + zsv(ji,jj) - zsv(ji,jj-1)  )
197               ! add it to the general tracer trends
198               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
199               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
200#if defined key_trdtra || defined key_trdmld
201               ! save the horizontal diffusive trends
202               ttrd(ji,jj,jk,3) = zta
203               strd(ji,jj,jk,3) = zsa
204#endif
205            END DO
206         END DO
207         !                                             ! ===============
208      END DO                                           ! Horizontal slab
209      !                                                ! ===============
210
211#if defined key_diaptr
212      ! "zonal" mean lateral diffusive heat and salt transport
213      == >> forced bug NOT implemented, ztv, zsv not 3D arrays
214      IF( MOD( kt, nf_ptr ) == 0 ) THEN
215# if defined key_s_coord || defined key_partial_steps
216         pht_ldf(:,:) = prt_vj( ztv(:,:,:) )
217         pst_ldf(:,:) = prt_vj( zsv(:,:,:) )
218# else
219         DO jk = 1, jpkm1
220            DO jj = 2, jpjm1
221               DO ji = fs_2, fs_jpim1   ! vector opt.
222                 ztv(ji,jj,jk) = ztv(ji,jj,jk) * fse3v(ji,jj,jk)
223                 zsv(ji,jj,jk) = zsv(ji,jj,jk) * fse3v(ji,jj,jk)
224               END DO
225            END DO
226         END DO
227         pht_ldf(:,:) = prt_vj( ztv(:,:,:) )
228         pst_ldf(:,:) = prt_vj( zsv(:,:,:) )
229# endif
230      ENDIF
231#endif
232
233   END SUBROUTINE tra_ldf_bilap
234
235   !!==============================================================================
236END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.