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

Last change on this file since 216 was 216, checked in by opalod, 19 years ago

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.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
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC tra_ldf_bilap   ! routine called by step.F90
27
28   !! * Substitutions
29#  include "domzgr_substitute.h90"
30#  include "ldftra_substitute.h90"
31#  include "ldfeiv_substitute.h90"
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !!   OPA 9.0 , LODYC-IPSL (2003)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38   
39   SUBROUTINE tra_ldf_bilap( kt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE tra_ldf_bilap  ***
42      !!
43      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
44      !!      trend and add it to the general trend of tracer equation.
45      !!
46      !! ** Method  :   4th order diffusive operator along model level surfaces
47      !!      evaluated using before fields (forward time scheme). The hor.
48      !!      diffusive trends of temperature (idem for salinity) is given by:
49      !!       * s-coordinate ('key_s_coord' defined), the vertical scale
50      !!      factors e3. are inside the derivatives:
51      !!      Laplacian of tb:
52      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
53      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
54      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
55      !!        zlt   = ahtt * zlt
56      !!        call to lbc_lnk
57      !!      Bilaplacian (laplacian of zlt):
58      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
59      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
60      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
61      !!      Laplacian of tb:
62      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(tb) ]
63      !!                              + dj-1[ e1v/e2v dj(tb) ] }
64      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
65      !!        zlt   = ahtt * zlt
66      !!        call to lbc_lnk
67      !!      Bilaplacian (laplacian of zlt):
68      !!         difft = 1/(e1t*e2t) {  di-1[ e2u/e1u di(zlt) ]
69      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
70      !!
71      !!      Add this trend to the general trend (ta,sa):
72      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
73      !!
74      !! ** Action : - Update (ta,sa) arrays with the before iso-level
75      !!               biharmonic mixing trend.
76      !!             - Save the trends in (ztdta,ztdsa) ('key_trdtra')
77      !!
78      !! History :
79      !!        !  91-11  (G. Madec)  Original code
80      !!        !  93-03  (M. Guyon)  symetrical conditions
81      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
82      !!        !  96-01  (G. Madec)  statement function for e3
83      !!        !  96-01  (M. Imbard)  mpp exchange
84      !!        !  97-07  (G. Madec)  optimization, and ahtt
85      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
86      !!   9.0  !  04-08  (C. Talandier) New trends organization
87      !!----------------------------------------------------------------------
88      !! * Modules used
89      USE oce           , ztu => ua,  &  ! use ua as workspace
90         &                ztv => va      ! use va as workspace
91
92      !! * Arguments
93      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
94
95      !! * Local declarations
96      INTEGER ::   ji, jj, jk             ! dummy loop indices
97#if defined key_partial_steps
98      INTEGER ::   iku, ikv               ! temporary integers
99#endif
100      REAL(wp) ::   zta, zsa              ! temporary scalars
101      REAL(wp), DIMENSION(jpi,jpj) ::   & 
102         zeeu, zeev, zbtr,              & ! workspace
103         zlt, zls
104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
105         zsu, zsv,                          & ! workspace arrays
106         ztdta, ztdsa
107      !!----------------------------------------------------------------------
108
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
112         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
113      ENDIF
114
115      ! Save ta and sa trends
116      IF( l_trdtra )   THEN
117         ztdta(:,:,:) = ta(:,:,:) 
118         ztdsa(:,:,:) = sa(:,:,:) 
119      ENDIF
120
121      !                                                ! ===============
122      DO jk = 1, jpkm1                                 ! Horizontal slab
123         !                                             ! ===============
124
125         ! 0. Initialization of metric arrays (for z- or s-coordinates)
126         ! ----------------------------------
127
128         DO jj = 1, jpjm1
129            DO ji = 1, fs_jpim1   ! vector opt.
130#if defined key_s_coord || defined key_partial_steps
131               ! s-coordinates, vertical scale factor are used
132               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
133               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
134               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
135#else
136               ! z-coordinates, no vertical scale factors
137               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
138               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
139               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
140#endif
141            END DO
142         END DO
143
144
145         ! 1. Laplacian
146         ! ------------
147
148         ! First derivative (gradient)
149         DO jj = 1, jpjm1
150            DO ji = 1, fs_jpim1   ! vector opt.
151               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
152               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
153               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
154               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
155            END DO
156         END DO
157#if defined key_partial_steps
158         DO jj = 1, jpj-1
159            DO ji = 1, jpi-1
160               ! last level
161               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
162               ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
163               IF( iku == jk ) THEN
164                  ztu(ji,jj,jk) = zeeu(ji,jj) * gtu(ji,jj)
165                  zsu(ji,jj,jk) = zeeu(ji,jj) * gsu(ji,jj)
166               ENDIF
167               IF( ikv == jk ) THEN
168                  ztv(ji,jj,jk) = zeev(ji,jj) * gtv(ji,jj)
169                  zsv(ji,jj,jk) = zeev(ji,jj) * gsv(ji,jj)
170               ENDIF
171            END DO
172         END DO
173#endif
174
175         ! Second derivative (divergence)
176         DO jj = 2, jpjm1
177            DO ji = fs_2, fs_jpim1   ! vector opt.
178               zlt(ji,jj) = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
179               zls(ji,jj) = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
180            END DO
181         END DO
182
183         ! Multiply by the eddy diffusivity coefficient
184         DO jj = 2, jpjm1
185            DO ji = fs_2, fs_jpim1   ! vector opt.
186               zlt(ji,jj) = fsahtt(ji,jj,jk) * zlt(ji,jj)
187               zls(ji,jj) = fsahtt(ji,jj,jk) * zls(ji,jj)
188            END DO
189         END DO
190
191         ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
192         CALL lbc_lnk( zlt, 'T', 1. )   ;    CALL lbc_lnk( zls, 'T', 1. )
193
194         ! 2. Bilaplacian
195         ! --------------
196
197         ! third derivative (gradient)
198         DO jj = 1, jpjm1
199            DO ji = 1, fs_jpim1   ! vector opt.
200               ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
201               zsu(ji,jj,jk) = zeeu(ji,jj) * ( zls(ji+1,jj  ) - zls(ji,jj) )
202               ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
203               zsv(ji,jj,jk) = zeev(ji,jj) * ( zls(ji  ,jj+1) - zls(ji,jj) )
204            END DO
205         END DO
206
207         ! fourth derivative (divergence) and add to the general tracer trend
208         DO jj = 2, jpjm1
209            DO ji = fs_2, fs_jpim1   ! vector opt.
210               ! horizontal diffusive trends
211               zta = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
212               zsa = zbtr(ji,jj) * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk) + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
213               ! add it to the general tracer trends
214               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
215               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
216            END DO
217         END DO
218         !                                             ! ===============
219      END DO                                           ! Horizontal slab
220      !                                                ! ===============
221
222      ! save the trends for diagnostic
223      ! save the horizontal diffusive trends
224      IF( l_trdtra )   THEN
225         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
226         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:)
227
228         CALL trd_mod(ztdta, ztdsa, jpttdldf, 'TRA', kt)
229      ENDIF
230
231      IF(l_ctl) THEN         ! print mean trends (used for debugging)
232         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
233         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
234         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
235         t_ctl = zta   ;   s_ctl = zsa
236      ENDIF
237
238      ! "zonal" mean lateral diffusive heat and salt transport
239      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
240# if defined key_s_coord || defined key_partial_steps
241         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
242         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
243# else
244         DO jk = 1, jpkm1
245            DO jj = 2, jpjm1
246               DO ji = fs_2, fs_jpim1   ! vector opt.
247                  ztv(ji,jj,jk) = ztv(ji,jj,jk) * fse3v(ji,jj,jk)
248                  zsv(ji,jj,jk) = zsv(ji,jj,jk) * fse3v(ji,jj,jk)
249               END DO
250            END DO
251         END DO
252         pht_ldf(:) = ptr_vj( ztv(:,:,:) )
253         pst_ldf(:) = ptr_vj( zsv(:,:,:) )
254# endif
255      ENDIF
256
257   END SUBROUTINE tra_ldf_bilap
258
259   !!==============================================================================
260END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.