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

Last change on this file since 1119 was 1119, checked in by cetlod, 16 years ago

style of all top namelist has been modified ; update modules to take it into account, see ticket:196

  • Property svn:executable set to *
File size: 10.2 KB
Line 
1MODULE trcldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  trcldf_bilap  ***
4   !! TOP :  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if defined key_top
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!----------------------------------------------------------------------
10   !!   trc_ldf_bilap : update the tracer trend with the horizontal diffusion
11   !!                   using a iso-level biharmonic operator
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce_trc         ! ocean dynamics and active tracers variables
15   USE trp_trc             ! ocean passive tracers variables
16   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
17   USE prtctl_trc      ! Print control for debbuging
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC trc_ldf_bilap   ! routine called by step.F90
24
25   !! * Substitutions
26#  include "top_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   TOP 1.0 , LOCEAN-IPSL (2005)
29   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcldf_bilap.F90,v 1.12 2007/10/12 09:26:30 opalod Exp $
30   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
31   !!----------------------------------------------------------------------
32
33CONTAINS
34   
35   SUBROUTINE trc_ldf_bilap( kt )
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE trc_ldf_bilap  ***
38      !!
39      !! ** Purpose :   Compute the before horizontal tracer tra diffusive
40      !!      trend and add it to the general trend of tracer equation.
41      !!
42      !! ** Method  :   4th order diffusive operator along model level surfaces
43      !!      evaluated using before fields (forward time scheme). The hor.
44      !!      diffusive trends of passive tracer is given by:
45      !!       * s-coordinate, the vertical scale
46      !!      factors e3. are inside the derivatives:
47      !!      Laplacian of trb:
48      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(trb) ]
49      !!                                  + dj-1[ e1v*e3v/e2v dj(trb) ]  }
50      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
51      !!        zlt   = ahtt * zlt
52      !!        call to lbc_lnk
53      !!      Bilaplacian (laplacian of zlt):
54      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
55      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
56      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
57      !!      Laplacian of trb:
58      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(trb) ]
59      !!                              + dj-1[ e1v/e2v dj(trb) ] }
60      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
61      !!        zlt   = ahtt * zlt
62      !!        call to lbc_lnk
63      !!      Bilaplacian (laplacian of zlt):
64      !!         difft = 1/(e1t*e2t) {  di-1[ e2u/e1u di(zlt) ]
65      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
66      !!
67      !!      Add this trend to the general trend tra :
68      !!         tra = tra + difft
69      !!
70      !! ** Action : - Update tra arrays with the before iso-level
71      !!               biharmonic mixing trend.
72      !!             - Save the trends in trtrd ('key_trc_diatrd')
73      !!
74      !! History :
75      !!        !  91-11  (G. Madec)  Original code
76      !!        !  93-03  (M. Guyon)  symetrical conditions
77      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
78      !!        !  96-01  (G. Madec)  statement function for e3
79      !!        !  96-01  (M. Imbard)  mpp exchange
80      !!        !  97-07  (G. Madec)  optimization, and ahtt
81      !!        !  00-05  (MA Foujols) add lbc for tracer trends
82      !!        !  00-10  (MA Foujols E. Kestenare) use passive tracer coefficient
83      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
84      !!   9.0  !  04-03  (C. Ethe )  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, jn             ! dummy loop indices
91      INTEGER ::   iku, ikv                   ! temporary integers
92
93      REAL(wp) ::   ztra     ! temporary scalars
94
95      REAL(wp), DIMENSION(jpi,jpj) ::   & 
96         zeeu, zeev, zbtr, zlt                 ! workspace
97      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
98         ztu, ztv                              ! workspace
99      CHARACTER (len=22) :: charout
100      !!----------------------------------------------------------------------
101
102      IF( kt == nittrc000 ) THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'trc_ldf_bilap : iso-level biharmonic operator'
105         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
106      ENDIF
107      !
108
109      DO jn = 1, jptra
110                                                          ! ===============
111         DO jk = 1, jpkm1                                 ! Horizontal slab
112            !                                             ! ===============
113
114            ! 0. Initialization of metric arrays (for z- or s-coordinates)
115            ! ----------------------------------
116
117            DO jj = 1, jpjm1
118               DO ji = 1, fs_jpim1   ! vector opt.
119#if ! defined key_zco
120                  ! s-coordinates, vertical scale factor are used
121                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
122                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
123                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
124#else
125                  ! z-coordinates, no vertical scale factors
126                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
127                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
128                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
129#endif
130               END DO
131            END DO
132
133
134            ! 1. Laplacian
135            ! ------------
136
137            ! First derivative (gradient)
138            DO jj = 1, jpjm1
139               DO ji = 1, fs_jpim1   ! vector opt.
140                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( trb(ji+1,jj  ,jk,jn) - trb(ji,jj,jk,jn) )
141                  ztv(ji,jj,jk) = zeev(ji,jj) * ( trb(ji  ,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
142               END DO
143            END DO
144
145            IF( ln_zps ) THEN
146               DO jj = 1, jpj-1
147                  DO ji = 1, jpi-1
148                     ! last level
149                     iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
150                     ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
151                     IF( iku == jk ) THEN
152                        ztu(ji,jj,jk) = zeeu(ji,jj) * gtru(ji,jj,jn)
153                     ENDIF
154                     IF( ikv == jk ) THEN
155                        ztv(ji,jj,jk) = zeev(ji,jj) * gtrv(ji,jj,jn)
156                     ENDIF
157                  END DO
158               END DO
159            ENDIF
160
161            ! Second derivative (divergence)
162            DO jj = 2, jpjm1
163               DO ji = fs_2, fs_jpim1   ! vector opt.
164                  zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
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) = fsahtrt(ji,jj,jk) * zlt(ji,jj)
172               END DO
173            END DO
174
175            ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
176            CALL lbc_lnk( zlt, '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,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
185                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
186               END DO
187            END DO
188
189            ! fourth derivative (divergence) and add to the general tracer trend
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1   ! vector opt.
192                  ! horizontal diffusive trends
193                  ztra = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
194                  ! add it to the general tracer trends
195                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
196#if defined key_trc_diatrd
197                  ! save the horizontal diffusive trends
198                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),4) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr(ji,jj)
199                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),5) = (  ztv(ji,jj,jk) - ztv(ji-1,jj,jk) ) * zbtr(ji,jj)
200#endif
201               END DO
202            END DO
203            !                                             ! ===============
204         END DO                                           ! Horizontal slab
205         !                                                ! ===============
206#if defined key_trc_diatrd
207         ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
208         IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),5), 'T', 1. ) 
209#endif
210      END DO
211
212     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
213         WRITE(charout, FMT="('ldf - bilap')")
214         CALL prt_ctl_trc_info(charout)
215         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
216      ENDIF
217
218   END SUBROUTINE trc_ldf_bilap
219
220#else
221   !!----------------------------------------------------------------------
222   !!   Default option                                         Empty module
223   !!----------------------------------------------------------------------
224CONTAINS
225   SUBROUTINE trc_ldf_bilap( kt ) 
226      INTEGER, INTENT(in) :: kt
227      WRITE(*,*) 'trc_ldf_bilap: You should not have seen this print! error?', kt
228   END SUBROUTINE trc_ldf_bilap
229#endif
230   !!==============================================================================
231END MODULE trcldf_bilap
Note: See TracBrowser for help on using the repository browser.