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

Last change on this file since 527 was 474, checked in by opalod, 18 years ago

nemo_v1_update_061: SM: end of ctl_stop + mpi optimization in _bilap
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 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   USE prtctl          ! Print control
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC tra_ldf_bilap   ! routine called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "ldftra_substitute.h90"
32#  include "ldfeiv_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LOCEAN-IPSL (2005)
36   !! $Header$
37   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
38   !!----------------------------------------------------------------------
39
40CONTAINS
41   
42   SUBROUTINE tra_ldf_bilap( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE tra_ldf_bilap  ***
45      !!
46      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
47      !!      trend and add it to the general trend of tracer equation.
48      !!
49      !! ** Method  :   4th order diffusive operator along model level surfaces
50      !!      evaluated using before fields (forward time scheme). The hor.
51      !!      diffusive trends of temperature (idem for salinity) is given by:
52      !!      Laplacian of tb:
53      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
54      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
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*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
60      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
61      !!      Note: if key_zco defined, e3t=e3u=e3v, they are simplified.
62      !!
63      !!      Add this trend to the general trend (ta,sa):
64      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
65      !!
66      !! ** Action : - Update (ta,sa) arrays with the before iso-level
67      !!               biharmonic mixing trend.
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      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
77      !!   9.0  !  04-08  (C. Talandier) New trends organization
78      !!        !  05-11  (G. Madec)  zps or sco as default option
79      !!----------------------------------------------------------------------
80      !! * Modules used
81      USE oce           , ztu => ua,  &  ! use ua as workspace
82         &                ztv => va      ! use va as workspace
83
84      !! * Arguments
85      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
86
87      !! * Local declarations
88      INTEGER ::   ji, jj, jk             ! dummy loop indices
89      INTEGER ::   iku, ikv               ! temporary integers
90      REAL(wp) ::   zta, zsa              ! temporary scalars
91      REAL(wp), DIMENSION(jpi,jpj) ::   & 
92         zeeu, zeev, zbtr                 ! 2D workspace arrays
93      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
94         zsu, zsv, zlt, zls               ! 3D workspace arrays
95      !!----------------------------------------------------------------------
96
97      IF( kt == nit000 ) THEN
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
100         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
101      ENDIF
102
103      !                                                ! ===============
104      DO jk = 1, jpkm1                                 ! Horizontal slab
105         !                                             ! ===============
106
107         ! 0. Initialization of metric arrays (for z- or s-coordinates)
108         ! ----------------------------------
109
110         IF( lk_zco ) THEN      ! z-coordinate (1D arrays): no vertical scale factors
111            DO jj = 1, jpjm1
112               DO ji = 1, fs_jpim1   ! vector opt.
113                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
114                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
115                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
116               END DO
117            END DO
118         ELSE                   ! All coordinates (3D arrays): vertical scale factor are used
119            DO jj = 1, jpjm1
120               DO ji = 1, fs_jpim1   ! vector opt.
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               END DO
125            END DO
126         ENDIF
127
128
129         ! 1. Laplacian
130         ! ------------
131
132         ! First derivative (gradient)
133         DO jj = 1, jpjm1
134            DO ji = 1, fs_jpim1   ! vector opt.
135               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
136               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
137               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
138               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
139            END DO
140         END DO
141         IF( ln_zps ) THEN      ! set gradient at partial step level
142            DO jj = 1, jpjm1
143               DO ji = 1, jpim1
144                  ! last level
145                  iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
146                  ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
147                  IF( iku == jk ) THEN
148                     ztu(ji,jj,jk) = zeeu(ji,jj) * gtu(ji,jj)
149                     zsu(ji,jj,jk) = zeeu(ji,jj) * gsu(ji,jj)
150                  ENDIF
151                  IF( ikv == jk ) THEN
152                     ztv(ji,jj,jk) = zeev(ji,jj) * gtv(ji,jj)
153                     zsv(ji,jj,jk) = zeev(ji,jj) * gsv(ji,jj)
154                  ENDIF
155               END DO
156            END DO
157         ENDIF
158
159         ! Second derivative (divergence)
160         DO jj = 2, jpjm1
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zlt(ji,jj,jk) = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
163               zls(ji,jj,jk) = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
164            END DO
165         END DO
166
167         ! Multiply by the eddy diffusivity coefficient
168         DO jj = 2, jpjm1
169            DO ji = fs_2, fs_jpim1   ! vector opt.
170               zlt(ji,jj,jk) = fsahtt(ji,jj,jk) * zlt(ji,jj,jk)
171               zls(ji,jj,jk) = fsahtt(ji,jj,jk) * zls(ji,jj,jk)
172            END DO
173         END DO
174      ENDDO
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      DO jk = 1, jpkm1
180     
181         ! 2. Bilaplacian
182         ! --------------
183
184         ! third derivative (gradient)
185         DO jj = 1, jpjm1
186            DO ji = 1, fs_jpim1   ! vector opt.
187               ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ,jk) - zlt(ji,jj,jk) )
188               zsu(ji,jj,jk) = zeeu(ji,jj) * ( zls(ji+1,jj  ,jk) - zls(ji,jj,jk) )
189               ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1,jk) - zlt(ji,jj,jk) )
190               zsv(ji,jj,jk) = zeev(ji,jj) * ( zls(ji  ,jj+1,jk) - zls(ji,jj,jk) )
191            END DO
192         END DO
193
194         ! fourth derivative (divergence) and add to the general tracer trend
195         DO jj = 2, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt.
197               ! horizontal diffusive trends
198               zta = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
199               zsa = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
200               ! add it to the general tracer trends
201               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
202               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
203            END DO
204         END DO
205         !                                             ! ===============
206      END DO                                           ! Horizontal slab
207      !                                                ! ===============
208
209      ! "zonal" mean lateral diffusive heat and salt transport
210      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
211         IF( lk_zco ) THEN      ! z-coordinate (1D arrays): multiply by the vertical scale factor
212            DO jk = 1, jpkm1
213               DO jj = 2, jpjm1
214                  DO ji = fs_2, fs_jpim1   ! vector opt.
215                     ztv(ji,jj,jk) = ztv(ji,jj,jk) * fse3v(ji,jj,jk)
216                     zsv(ji,jj,jk) = zsv(ji,jj,jk) * fse3v(ji,jj,jk)
217                  END DO
218               END DO
219            END DO
220         ENDIF
221         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
222         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
223      ENDIF
224
225   END SUBROUTINE tra_ldf_bilap
226
227   !!==============================================================================
228END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.