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 tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/traldf_bilap.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
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   !! $Id: traldf_bilap.F90 1152 2008-06-26 14:11:13Z rblod $
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
93         zlt, zls
94      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
95         zsu, zsv                         ! 3D workspace
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
102      ENDIF
103
104
105      !                                                ! ===============
106      DO jk = 1, jpkm1                                 ! Horizontal slab
107         !                                             ! ===============
108
109         ! 0. Initialization of metric arrays (for z- or s-coordinates)
110         ! ----------------------------------
111
112         IF( lk_zco ) THEN      ! z-coordinate (1D arrays): no vertical scale factors
113            DO jj = 1, jpjm1
114               DO ji = 1, fs_jpim1   ! vector opt.
115                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
116                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
117                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
118               END DO
119            END DO
120         ELSE                   ! All coordinates (3D arrays): vertical scale factor are used
121            DO jj = 1, jpjm1
122               DO ji = 1, fs_jpim1   ! vector opt.
123                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
124                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
125                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
126               END DO
127            END DO
128         ENDIF
129
130
131         ! 1. Laplacian
132         ! ------------
133
134         ! First derivative (gradient)
135         DO jj = 1, jpjm1
136            DO ji = 1, fs_jpim1   ! vector opt.
137               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
138               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
139               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
140               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
141            END DO
142         END DO
143         IF( ln_zps ) THEN      ! set gradient at partial step level
144            DO jj = 1, jpjm1
145               DO ji = 1, jpim1
146                  ! last level
147                  iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
148                  ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
149                  IF( iku == jk ) THEN
150                     ztu(ji,jj,jk) = zeeu(ji,jj) * gtu(ji,jj)
151                     zsu(ji,jj,jk) = zeeu(ji,jj) * gsu(ji,jj)
152                  ENDIF
153                  IF( ikv == jk ) THEN
154                     ztv(ji,jj,jk) = zeev(ji,jj) * gtv(ji,jj)
155                     zsv(ji,jj,jk) = zeev(ji,jj) * gsv(ji,jj)
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               zls(ji,jj) = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
166            END DO
167         END DO
168
169         ! Multiply by the eddy diffusivity coefficient
170         DO jj = 2, jpjm1
171            DO ji = fs_2, fs_jpim1   ! vector opt.
172               zlt(ji,jj) = fsahtt(ji,jj,jk) * zlt(ji,jj)
173               zls(ji,jj) = fsahtt(ji,jj,jk) * zls(ji,jj)
174            END DO
175         END DO
176
177         ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
178         CALL lbc_lnk( zlt, 'T', 1. )   ;    CALL lbc_lnk( zls, 'T', 1. )
179
180         ! 2. Bilaplacian
181         ! --------------
182
183         ! third derivative (gradient)
184         DO jj = 1, jpjm1
185            DO ji = 1, fs_jpim1   ! vector opt.
186               ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
187               zsu(ji,jj,jk) = zeeu(ji,jj) * ( zls(ji+1,jj  ) - zls(ji,jj) )
188               ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
189               zsv(ji,jj,jk) = zeev(ji,jj) * ( zls(ji  ,jj+1) - zls(ji,jj) )
190            END DO
191         END DO
192
193         ! fourth derivative (divergence) and add to the general tracer trend
194         DO jj = 2, jpjm1
195            DO ji = fs_2, fs_jpim1   ! vector opt.
196               ! horizontal diffusive trends
197               zta = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
198               zsa = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
199               ! add it to the general tracer trends
200               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
201               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
202            END DO
203         END DO
204         !                                             ! ===============
205      END DO                                           ! Horizontal slab
206      !                                                ! ===============
207
208      ! "zonal" mean lateral diffusive heat and salt transport
209      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
210         IF( lk_zco ) THEN      ! z-coordinate (1D arrays): multiply by the vertical scale factor
211            DO jk = 1, jpkm1
212               DO jj = 2, jpjm1
213                  DO ji = fs_2, fs_jpim1   ! vector opt.
214                     ztv(ji,jj,jk) = ztv(ji,jj,jk) * fse3v(ji,jj,jk)
215                     zsv(ji,jj,jk) = zsv(ji,jj,jk) * fse3v(ji,jj,jk)
216                  END DO
217               END DO
218            END DO
219         ENDIF
220         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
221         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
222      ENDIF
223
224   END SUBROUTINE tra_ldf_bilap
225
226   !!==============================================================================
227END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.