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.
trcldf_bilap.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcldf_bilap.F90 @ 197

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

CL + CE : NEMO TRC_SRC start
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.1 KB
Line 
1MODULE trcldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  trcldf_bilap  ***
4   !! Ocean passive tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_ldf_bilap : update the tracer trend with the horizontal diffusion
9   !!                   using a iso-level biharmonic operator
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce_trc         ! ocean dynamics and active tracers
13   USE trc             ! ocean space and time domain
14   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
15
16   IMPLICIT NONE
17   PRIVATE
18
19   !! * Routine accessibility
20   PUBLIC trc_ldf_bilap   ! routine called by step.F90
21
22   !! * Substitutions
23#  include "passivetrc_substitute.h90"
24   !!----------------------------------------------------------------------
25   !!   OPA 9.0 , LODYC-IPSL (2003)
26   !!----------------------------------------------------------------------
27
28CONTAINS
29   
30   SUBROUTINE trc_ldf_bilap( kt )
31      !!----------------------------------------------------------------------
32      !!                  ***  ROUTINE trc_ldf_bilap  ***
33      !!
34      !! ** Purpose :   Compute the before horizontal tracer tra diffusive
35      !!      trend and add it to the general trend of tracer equation.
36      !!
37      !! ** Method  :   4th order diffusive operator along model level surfaces
38      !!      evaluated using before fields (forward time scheme). The hor.
39      !!      diffusive trends of passive tracer is given by:
40      !!       * s-coordinate ('key_s_coord' defined), the vertical scale
41      !!      factors e3. are inside the derivatives:
42      !!      Laplacian of trb:
43      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(trb) ]
44      !!                                  + dj-1[ e1v*e3v/e2v dj(trb) ]  }
45      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
46      !!        zlt   = ahtt * zlt
47      !!        call to lbc_lnk
48      !!      Bilaplacian (laplacian of zlt):
49      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
50      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
51      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
52      !!      Laplacian of trb:
53      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(trb) ]
54      !!                              + dj-1[ e1v/e2v dj(trb) ] }
55      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
56      !!        zlt   = ahtt * zlt
57      !!        call to lbc_lnk
58      !!      Bilaplacian (laplacian of zlt):
59      !!         difft = 1/(e1t*e2t) {  di-1[ e2u/e1u di(zlt) ]
60      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
61      !!
62      !!      Add this trend to the general trend tra :
63      !!         tra = tra + difft
64      !!
65      !! ** Action : - Update tra arrays with the before iso-level
66      !!               biharmonic mixing trend.
67      !!             - Save the trends in trtrd ('key_trc_diatrd')
68      !!
69      !! History :
70      !!        !  91-11  (G. Madec)  Original code
71      !!        !  93-03  (M. Guyon)  symetrical conditions
72      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
73      !!        !  96-01  (G. Madec)  statement function for e3
74      !!        !  96-01  (M. Imbard)  mpp exchange
75      !!        !  97-07  (G. Madec)  optimization, and ahtt
76      !!        !  00-05  (MA Foujols) add lbc for tracer trends
77      !!        !  00-10  (MA Foujols E. Kestenare) use passive tracer coefficient
78      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
79      !!   9.0  !  04-03  (C. Ethe )  F90: Free form and module
80      !!----------------------------------------------------------------------
81      !! * Arguments
82      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
83
84      !! * Local declarations
85      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
86
87      REAL(wp) ::   ztra     ! temporary scalars
88
89      REAL(wp), DIMENSION(jpi,jpj) ::   & 
90         zeeu, zeev, zbtr,             &  ! workspace
91         zlt, ztu, ztv
92      !!----------------------------------------------------------------------
93
94      IF( kt == nittrc000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'trc_ldf_bilap : iso-level biharmonic operator'
97         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
98      ENDIF
99      !
100
101      DO jn = 1, jptra
102                                                          ! ===============
103         DO jk = 1, jpkm1                                 ! Horizontal slab
104            !                                             ! ===============
105
106            ! 0. Initialization of metric arrays (for z- or s-coordinates)
107            ! ----------------------------------
108
109            DO jj = 1, jpjm1
110               DO ji = 1, fs_jpim1   ! vector opt.
111#if defined key_s_coord || defined key_partial_steps
112                  ! s-coordinates, vertical scale factor are used
113                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
114                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
115                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
116#else
117                  ! z-coordinates, no vertical scale factors
118                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
119                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
120                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
121#endif
122               END DO
123            END DO
124
125
126            ! 1. Laplacian
127            ! ------------
128
129            ! First derivative (gradient)
130            DO jj = 1, jpjm1
131               DO ji = 1, fs_jpim1   ! vector opt.
132                  ztu(ji,jj) = zeeu(ji,jj) * ( trb(ji+1,jj  ,jk,jn) - trb(ji,jj,jk,jn) )
133                  ztv(ji,jj) = zeev(ji,jj) * ( trb(ji  ,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
134               END DO
135            END DO
136
137
138            ! Second derivative (divergence)
139            DO jj = 2, jpjm1
140               DO ji = fs_2, fs_jpim1   ! vector opt.
141                  zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj) - ztu(ji-1,jj) + ztv(ji,jj) - ztv(ji,jj-1)  )
142               END DO
143            END DO
144
145            ! Multiply by the eddy diffusivity coefficient
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  zlt(ji,jj) = fsahtrt(ji,jj,jk) * zlt(ji,jj)
149               END DO
150            END DO
151
152            ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
153            CALL lbc_lnk( zlt, 'T', 1. ) 
154
155            ! 2. Bilaplacian
156            ! --------------
157
158            ! third derivative (gradient)
159            DO jj = 1, jpjm1
160               DO ji = 1, fs_jpim1   ! vector opt.
161                  ztu(ji,jj) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
162                  ztv(ji,jj) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
163               END DO
164            END DO
165
166            ! fourth derivative (divergence) and add to the general tracer trend
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.
169                  ! horizontal diffusive trends
170                  ztra = zbtr(ji,jj) * (  ztu(ji,jj) - ztu(ji-1,jj) + ztv(ji,jj) - ztv(ji,jj-1)  )
171                  ! add it to the general tracer trends
172                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
173#if defined key_trc_diatrd
174                  ! save the horizontal diffusive trends
175                  trtrd(ji,jj,jk,jn,4) = (  ztu(ji,jj) - ztu(ji-1,jj) ) * zbtr(ji,jj)
176                  trtrd(ji,jj,jk,jn,5) = (  ztv(ji,jj) - ztv(ji-1,jj) ) * zbtr(ji,jj)
177#endif
178               END DO
179            END DO
180            !                                             ! ===============
181         END DO                                           ! Horizontal slab
182         !                                                ! ===============
183#if defined key_trc_diatrd
184         ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
185         CALL lbc_lnk( trtrd(1,1,1,jn,5), 'T', 1. ) 
186#endif
187         IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
188            ztra = SUM( tra(2:jpim1,2:jpjm1,1:jpkm1,jn) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
189            WRITE(numout,*) ' trc/ldf  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn)
190            tra_ctl(jn) = ztra 
191         ENDIF
192
193      END DO
194
195   END SUBROUTINE trc_ldf_bilap
196
197#else
198   !!----------------------------------------------------------------------
199   !!   Default option                                         Empty module
200   !!----------------------------------------------------------------------
201CONTAINS
202   SUBROUTINE trc_ldf_bilap( kt ) 
203      INTEGER, INTENT(in) :: kt
204      WRITE(*,*) 'trc_ldf_bilap: You should not have seen this print! error?', kt
205   END SUBROUTINE trc_ldf_bilap
206#endif
207   !!==============================================================================
208END MODULE trcldf_bilap
Note: See TracBrowser for help on using the repository browser.