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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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