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

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

CT : UPDATE082 : Finalization of the poleward transport diagnostics
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.3 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 diaptr          ! 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      !! * Modules used
87      USE oce           , ztu => ua,  &  ! use ua as workspace
88         &                ztv => va      ! use va as workspace
89
90      !! * Arguments
91      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
92
93      !! * Local declarations
94      INTEGER ::   ji, jj, jk             ! dummy loop indices
95#if defined key_partial_steps
96      INTEGER ::   iku, ikv               ! temporary integers
97#endif
98      REAL(wp) ::   zta, zsa    ! temporary scalars
99      REAL(wp), DIMENSION(jpi,jpj) ::   & 
100         zeeu, zeev, zbtr,             &  ! workspace
101         zlt, zls
102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
103         zsu, zsv
104      !!----------------------------------------------------------------------
105
106      IF( kt == nit000 ) THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
109         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
110      ENDIF
111      !                                                ! ===============
112      DO jk = 1, jpkm1                                 ! Horizontal slab
113         !                                             ! ===============
114
115         ! 0. Initialization of metric arrays (for z- or s-coordinates)
116         ! ----------------------------------
117
118         DO jj = 1, jpjm1
119            DO ji = 1, fs_jpim1   ! vector opt.
120#if defined key_s_coord || defined key_partial_steps
121               ! s-coordinates, vertical scale factor are used
122               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
123               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
124               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
125#else
126               ! z-coordinates, no vertical scale factors
127               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
128               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
129               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
130#endif
131            END DO
132         END DO
133
134
135         ! 1. Laplacian
136         ! ------------
137
138         ! First derivative (gradient)
139         DO jj = 1, jpjm1
140            DO ji = 1, fs_jpim1   ! vector opt.
141               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
142               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
143               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
144               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
145            END DO
146         END DO
147#if defined key_partial_steps
148         DO jj = 1, jpj-1
149            DO ji = 1, jpi-1
150               ! last level
151               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
152               ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
153               IF( iku == jk ) THEN
154                  ztu(ji,jj,jk) = zeeu(ji,jj) * gtu(ji,jj)
155                  zsu(ji,jj,jk) = zeeu(ji,jj) * gsu(ji,jj)
156               ENDIF
157               IF( ikv == jk ) THEN
158                  ztv(ji,jj,jk) = zeev(ji,jj) * gtv(ji,jj)
159                  zsv(ji,jj,jk) = zeev(ji,jj) * gsv(ji,jj)
160               ENDIF
161            END DO
162         END DO
163#endif
164
165         ! Second derivative (divergence)
166         DO jj = 2, jpjm1
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168               zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
169               zls(ji,jj) = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
170            END DO
171         END DO
172
173         ! Multiply by the eddy diffusivity coefficient
174         DO jj = 2, jpjm1
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               zlt(ji,jj) = fsahtt(ji,jj,jk) * zlt(ji,jj)
177               zls(ji,jj) = fsahtt(ji,jj,jk) * zls(ji,jj)
178            END DO
179         END DO
180
181         ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
182         CALL lbc_lnk( zlt, 'T', 1. )   ;    CALL lbc_lnk( zls, 'T', 1. )
183
184         ! 2. Bilaplacian
185         ! --------------
186
187         ! third derivative (gradient)
188         DO jj = 1, jpjm1
189            DO ji = 1, fs_jpim1   ! vector opt.
190               ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
191               zsu(ji,jj,jk) = zeeu(ji,jj) * ( zls(ji+1,jj  ) - zls(ji,jj) )
192               ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
193               zsv(ji,jj,jk) = zeev(ji,jj) * ( zls(ji  ,jj+1) - zls(ji,jj) )
194            END DO
195         END DO
196
197         ! fourth derivative (divergence) and add to the general tracer trend
198         DO jj = 2, jpjm1
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               ! horizontal diffusive trends
201               zta = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
202               zsa = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
203               ! add it to the general tracer trends
204               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
205               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
206#if defined key_trdtra || defined key_trdmld
207               ! save the horizontal diffusive trends
208               ttrd(ji,jj,jk,3) = zta
209               strd(ji,jj,jk,3) = zsa
210#endif
211            END DO
212         END DO
213         !                                             ! ===============
214      END DO                                           ! Horizontal slab
215      !                                                ! ===============
216
217      ! "zonal" mean lateral diffusive heat and salt transport
218      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
219# if defined key_s_coord || defined key_partial_steps
220         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
221         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
222# else
223         DO jk = 1, jpkm1
224            DO jj = 2, jpjm1
225               DO ji = fs_2, fs_jpim1   ! vector opt.
226                  ztv(ji,jj,jk) = ztv(ji,jj,jk) * fse3v(ji,jj,jk)
227                  zsv(ji,jj,jk) = zsv(ji,jj,jk) * fse3v(ji,jj,jk)
228               END DO
229            END DO
230         END DO
231         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
232         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
233# endif
234      ENDIF
235
236   END SUBROUTINE tra_ldf_bilap
237
238   !!==============================================================================
239END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.