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 @ 74

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

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