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.
dynldf_bilap_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynldf_bilap_tam.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007
  • Property svn:executable set to *
File size: 19.9 KB
Line 
1MODULE dynldf_bilap_tam
2#ifdef key_tam
3   !!===========================================================================
4   !!                       ***  MODULE  dynldf_bilap_tam  ***
5   !! Ocean dynamics:  lateral viscosity trend
6   !!                  Tangent and Adjoint Module
7   !!===========================================================================
8
9   !!---------------------------------------------------------------------------
10   !!   dyn_ldf_bilap_tan : update the momentum trend with the lateral diffusion
11   !!                       using an iso-level bilaplacian operator (tangent)
12   !!   dyn_ldf_bilap_adj : update the momentum trend with the lateral diffusion
13   !!                       using an iso-level bilaplacian operator (adjoint)
14   !!---------------------------------------------------------------------------
15   !! * Modules used
16   USE lbclnk
17   USE lbclnk_tam
18   USE par_oce
19   USE oce_tam
20   USE ldfdyn_oce
21   USE dom_oce
22   USE in_out_manager
23   USE wrk_nemo        ! Memory Allocation
24   USE timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC dyn_ldf_bilap_tan  ! called by dynldf_tam.F90
31   PUBLIC dyn_ldf_bilap_adj  ! called by dynldf_tam.F90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "ldfdyn_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE dyn_ldf_bilap_tan( kt )
42      !!----------------------------------------------------------------------
43      !!                     ***  ROUTINE dyn_ldf_bilap_tan  ***
44      !!
45      !! ** Purpose :   Compute the before trend of the lateral momentum
46      !!      diffusion and add it to the general trend of momentum equation.
47      !!
48      !! ** Method  :   The before horizontal momentum diffusion trend is a
49      !!      bi-harmonic operator (bilaplacian type) which separates the
50      !!      divergent and rotational parts of the flow.
51      !!      Its horizontal components are computed as follow:
52      !!      laplacian:
53      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
54      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
55      !!      third derivative:
56      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
57      !!          zlu = ahmu * zlu
58      !!          zlv = ahmv * zlv
59      !!       * curl and divergence of the laplacian
60      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
61      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
62      !!      bilaplacian:
63      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
64      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
65      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
66      !!      rotational part of the diffusion are simplified
67      !!      Add this before trend to the general trend (ua,va):
68      !!            (ua,va) = (ua,va) + (diffu,diffv)
69      !!      'key_trddyn' defined: the two components of the horizontal
70      !!                               diffusion trend are saved.
71      !!
72      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
73      !!               mixing trend.
74      !!
75      !! History :
76      !!        !  90-09  (G. Madec)  Original code
77      !!        !  91-11  (G. Madec)
78      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon)
79      !!        !  96-01  (G. Madec)  statement function for e3
80      !!        !  97-07  (G. Madec)  lbc calls
81      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
82      !!   9.0  !  04-08  (C. Talandier) New trends organization
83      !! History of the tangent routine
84      !!   9.0  !  09-12 (F. Vigilant) tangent of 9.0
85      !!   3.4  !  11-07 (P.-A. Bouttier) 3.4 version
86      !!----------------------------------------------------------------------
87      !! * Arguments
88      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
89      !! * Local declarations
90      INTEGER  ::   ji, jj, jk                ! dummy loop indices
91      REAL(wp) ::   &
92         zuatl, zvatl, zbt, ze2u, ze2v        ! temporary scalars
93      REAL(wp), POINTER, DIMENSION(:,:) ::       &
94         zcutl, zcvtl                         ! temporary workspace
95      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
96         zuftl, zuttl, zlutl, zlvtl           ! temporary workspace
97      !!----------------------------------------------------------------------
98      !
99      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap_tan')
100      !
101      CALL wrk_alloc( jpi, jpj,      zcutl, zcvtl           )
102      CALL wrk_alloc( jpi, jpj, jpk, zuftl, zuttl, zlutl, zlvtl )
103      !
104      IF( kt == nit000 ) THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap_tan: iso-level bilaplacien operator'
107         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ '
108      ENDIF
109
110      zuftl(:,:,:) = 0.0_wp
111      zuttl(:,:,:) = 0.0_wp
112      zlutl(:,:,:) = 0.0_wp
113      zlvtl(:,:,:) = 0.0_wp
114      !                                                ! ===============
115      DO jk = 1, jpkm1                                 ! Horizontal slab
116         !                                             ! ===============
117         ! Laplacian
118         ! ---------
119
120         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
121            zuftl(:,:,jk) = rotb_tl(:,:,jk) * fse3f(:,:,jk)
122            DO jj = 2, jpjm1
123               DO ji = fs_2, fs_jpim1   ! vector opt.
124                  zlutl(ji,jj,jk) = - ( zuftl(ji,jj,jk) - zuftl(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
125                     &         + ( hdivb_tl(ji+1,jj,jk) - hdivb_tl(ji,jj,jk) ) / e1u(ji,jj)
126
127                  zlvtl(ji,jj,jk) = + ( zuftl(ji,jj,jk) - zuftl(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
128                     &         + ( hdivb_tl(ji,jj+1,jk) - hdivb_tl(ji,jj,jk) ) / e2v(ji,jj)
129               END DO
130            END DO
131         ELSE                            ! z-coordinate - full step
132            DO jj = 2, jpjm1
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134                  zlutl(ji,jj,jk) = - ( rotb_tl (ji  ,jj,jk) - rotb_tl (ji,jj-1,jk) ) / e2u(ji,jj)   &
135                     &         + ( hdivb_tl(ji+1,jj,jk) - hdivb_tl(ji,jj  ,jk) ) / e1u(ji,jj)
136
137                  zlvtl(ji,jj,jk) = + ( rotb_tl (ji,jj  ,jk) - rotb_tl (ji-1,jj,jk) ) / e1v(ji,jj)   &
138                     &         + ( hdivb_tl(ji,jj+1,jk) - hdivb_tl(ji  ,jj,jk) ) / e2v(ji,jj)
139               END DO
140            END DO
141         ENDIF
142      ENDDO
143
144      ! Boundary conditions on the laplacian  (zlu,zlv)
145      CALL lbc_lnk( zlutl, 'U', -1.0_wp )
146      CALL lbc_lnk( zlvtl, 'V', -1.0_wp )
147
148      DO jk = 1, jpkm1
149
150         ! Third derivative
151         ! ----------------
152
153         ! Multiply by the eddy viscosity coef. (at u- and v-points)
154         zlutl(:,:,jk) = zlutl(:,:,jk) * fsahmu(:,:,jk)
155         zlvtl(:,:,jk) = zlvtl(:,:,jk) * fsahmv(:,:,jk)
156
157         ! Contravariant "laplacian"
158         zcutl(:,:) = e1u(:,:) * zlutl(:,:,jk)
159         zcvtl(:,:) = e2v(:,:) * zlvtl(:,:,jk)
160
161         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
162         DO jj = 1, jpjm1
163            DO ji = 1, fs_jpim1   ! vector opt.
164               zuftl(ji,jj,jk) = fmask(ji,jj,jk) * (  zcvtl(ji+1,jj  ) - zcvtl(ji,jj)      &
165                  &                            - zcutl(ji  ,jj+1) + zcutl(ji,jj)  )   &
166
167                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
168            END DO
169         END DO
170
171         ! Laplacian Horizontal fluxes
172         DO jj = 1, jpjm1
173            DO ji = 1, fs_jpim1   ! vector opt.
174               zlutl(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlutl(ji,jj,jk)
175               zlvtl(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlvtl(ji,jj,jk)
176            END DO
177         END DO
178
179         ! Laplacian divergence
180         DO jj = 2, jpj
181            DO ji = fs_2, jpi   ! vector opt.
182               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
183               zuttl(ji,jj,jk) = (  zlutl(ji,jj,jk) - zlutl(ji-1,jj  ,jk)   &
184                  &             + zlvtl(ji,jj,jk) - zlvtl(ji  ,jj-1,jk) ) / zbt
185            END DO
186         END DO
187      END DO
188
189      ! boundary conditions on the laplacian curl and div (zuf,zut)
190!!bug gm no need to do this 2 following lbc...
191      CALL lbc_lnk( zuftl, 'F', 1.0_wp )
192      CALL lbc_lnk( zuttl, 'T', 1.0_wp )
193
194      DO jk = 1, jpkm1
195
196         ! Bilaplacian
197         ! -----------
198
199         DO jj = 2, jpjm1
200            DO ji = fs_2, fs_jpim1   ! vector opt.
201               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
202               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
203               ! horizontal biharmonic diffusive trends
204               zuatl = - ( zuftl(ji  ,jj,jk) - zuftl(ji,jj-1,jk) ) / ze2u   &
205                  &  + ( zuttl(ji+1,jj,jk) - zuttl(ji,jj  ,jk) ) / e1u(ji,jj)
206               zvatl = + ( zuftl(ji,jj  ,jk) - zuftl(ji-1,jj,jk) ) / ze2v   &
207                  &  + ( zuttl(ji,jj+1,jk) - zuttl(ji  ,jj,jk) ) / e2v(ji,jj)
208               ! add it to the general momentum trends
209               ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zuatl
210               va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zvatl
211            END DO
212         END DO
213
214         !                                             ! ===============
215      END DO                                           !   End of slab
216      !                                                ! ===============
217      CALL wrk_dealloc( jpi, jpj,      zcutl, zcvtl               )
218      CALL wrk_dealloc( jpi, jpj, jpk, zuftl, zuttl, zlutl, zlvtl )
219      !
220      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap_tan')
221      !
222   END SUBROUTINE dyn_ldf_bilap_tan
223
224
225   SUBROUTINE dyn_ldf_bilap_adj( kt )
226      !!----------------------------------------------------------------------
227      !!                     ***  ROUTINE dyn_ldf_bilap_adj  ***
228      !!
229      !! ** Purpose :   Compute the before trend of the lateral momentum
230      !!      diffusion and add it to the general trend of momentum equation.
231      !!
232      !! ** Method  :   The before horizontal momentum diffusion trend is a
233      !!      bi-harmonic operator (bilaplacian type) which separates the
234      !!      divergent and rotational parts of the flow.
235      !!      Its horizontal components are computed as follow:
236      !!      laplacian:
237      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
238      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
239      !!      third derivative:
240      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
241      !!          zlu = ahmu * zlu
242      !!          zlv = ahmv * zlv
243      !!       * curl and divergence of the laplacian
244      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
245      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
246      !!      bilaplacian:
247      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
248      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
249      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
250      !!      rotational part of the diffusion are simplified
251      !!      Add this before trend to the general trend (ua,va):
252      !!            (ua,va) = (ua,va) + (diffu,diffv)
253      !!      'key_trddyn' defined: the two components of the horizontal
254      !!                               diffusion trend are saved.
255      !!
256      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
257      !!               mixing trend.
258      !!
259      !! History :
260      !!        !  90-09  (G. Madec)  Original code
261      !!        !  91-11  (G. Madec)
262      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon)
263      !!        !  96-01  (G. Madec)  statement function for e3
264      !!        !  97-07  (G. Madec)  lbc calls
265      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
266      !!   9.0  !  04-08  (C. Talandier) New trends organization
267      !! History of the adjoint routine
268      !!   3.4  !  12-07  (P.-A. Bouttier) 3.4 version
269      !!----------------------------------------------------------------------
270      !! * Arguments
271      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
272      !! * Local declarations
273      INTEGER  ::   ji, jj, jk            ! dummy loop indices
274      REAL(wp) ::   &
275         zuaad, zvaad, zbt, ze2u, ze2v        ! temporary scalars
276      REAL(wp), POINTER, DIMENSION(:,:) ::       &
277         zcuad, zcvad                         ! temporary workspace
278      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
279         zufad, zutad, zluad, zlvad           ! temporary workspace
280      !!----------------------------------------------------------------------
281      !
282      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap_adj')
283      !
284      CALL wrk_alloc( jpi, jpj,      zcuad, zcvad           )
285      CALL wrk_alloc( jpi, jpj, jpk, zufad, zutad, zluad, zlvad )
286      !
287      IF( kt == nitend ) THEN
288         IF(lwp) WRITE(numout,*)
289         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap_adj: bilaplacien operator'
290         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ '
291      ENDIF
292
293      zuaad = 0.0_wp
294      zvaad = 0.0_wp
295      zufad(:,:,:) = 0.0_wp
296      zutad(:,:,:) = 0.0_wp
297      zluad(:,:,:) = 0.0_wp
298      zlvad(:,:,:) = 0.0_wp
299      zcvad(:,:) = 0.0_wp
300      zcuad(:,:) = 0.0_wp
301
302      DO jk = 1, jpkm1
303         ! Bilaplacian
304         ! -----------
305         DO jj = jpjm1, 2, -1
306            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
307               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
308               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
309               ! add it to the general momentum trends
310               zvaad = zvaad + va_ad(ji,jj,jk)
311               zuaad = zuaad + ua_ad(ji,jj,jk)
312               ! horizontal biharmonic diffusive trends
313               zufad(ji  ,jj  ,jk) = zufad(ji  ,jj  ,jk) + zvaad / ze2v
314               zufad(ji-1,jj  ,jk) = zufad(ji-1,jj  ,jk) - zvaad / ze2v
315               zutad(ji  ,jj  ,jk) = zutad(ji  ,jj  ,jk) - zvaad / e2v(ji,jj)
316               zutad(ji  ,jj+1,jk) = zutad(ji  ,jj+1,jk) + zvaad / e2v(ji,jj)
317               !
318               zufad(ji  ,jj  ,jk) = zufad(ji  ,jj  ,jk) - zuaad / ze2u
319               zufad(ji  ,jj-1,jk) = zufad(ji  ,jj-1,jk) + zuaad / ze2u
320               zutad(ji  ,jj  ,jk) = zutad(ji  ,jj  ,jk) - zuaad / e1u(ji,jj)
321               zutad(ji+1,jj  ,jk) = zutad(ji+1,jj  ,jk) + zuaad / e1u(ji,jj)
322               !
323               zuaad = 0.0_wp
324               zvaad = 0.0_wp
325            END DO
326         END DO
327         !                                             ! ===============
328      END DO                                           !   End of slab
329      !                                                ! ===============
330      ! boundary conditions on the laplacian curl and div (zuf,zut)
331!!bug gm no need to do this 2 following lbc...
332      CALL lbc_lnk_adj( zutad, 'T', 1.0_wp )
333      CALL lbc_lnk_adj( zufad, 'F', 1.0_wp )
334      DO jk = 1, jpkm1
335         ! Third derivative
336         ! ----------------
337         ! Laplacian divergence
338         DO jj = jpj, 2, -1
339            DO ji = jpi, fs_2, -1   ! vector opt.
340               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
341               zluad(ji  ,jj  ,jk) = zluad(ji  ,jj  ,jk) + zutad(ji,jj,jk) / zbt
342               zluad(ji-1,jj  ,jk) = zluad(ji-1,jj  ,jk) - zutad(ji,jj,jk) / zbt
343               zlvad(ji  ,jj  ,jk) = zlvad(ji  ,jj  ,jk) + zutad(ji,jj,jk) / zbt
344               zlvad(ji  ,jj-1,jk) = zlvad(ji  ,jj-1,jk) - zutad(ji,jj,jk) / zbt
345               zutad(ji,jj,jk) = 0.0_wp
346            END DO
347         END DO
348         ! Laplacian Horizontal fluxes
349         DO jj = jpjm1, 1, -1
350            DO ji = fs_jpim1, 1, -1   ! vector opt.
351               zluad(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zluad(ji,jj,jk)
352               zlvad(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlvad(ji,jj,jk)
353            END DO
354         END DO
355         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
356         DO jj = jpjm1, 1, -1
357            DO ji = fs_jpim1, 1, -1   ! vector opt.
358               zufad(ji,jj,jk) = fmask(ji,jj,jk) * zufad(ji,jj,jk) * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
359               zcvad(ji  ,jj  ) = zcvad(ji  ,jj  ) - zufad(ji,jj,jk)
360               zcvad(ji+1,jj  ) = zcvad(ji+1,jj  ) + zufad(ji,jj,jk)
361               zcuad(ji  ,jj  ) = zcuad(ji  ,jj  ) + zufad(ji,jj,jk)
362               zcuad(ji  ,jj+1) = zcuad(ji  ,jj+1) - zufad(ji,jj,jk)
363               zufad(ji,jj,jk) = 0.0_wp
364            END DO
365         END DO
366         ! Contravariant "laplacian"
367         DO jj = 1, jpj
368            DO ji = 1, jpi
369               zlvad(ji,jj,jk) = zlvad(ji,jj,jk) + e2v(ji,jj) * zcvad(ji,jj)
370               zluad(ji,jj,jk) = zluad(ji,jj,jk) + e1u(ji,jj) * zcuad(ji,jj)
371               zcvad(ji,jj) = 0.0_wp
372               zcuad(ji,jj) = 0.0_wp
373            END DO
374         END DO
375         ! Multiply by the eddy viscosity coef. (at u- and v-points)
376         zluad(:,:,jk) = zluad(:,:,jk) * fsahmu(:,:,jk)
377         zlvad(:,:,jk) = zlvad(:,:,jk) * fsahmv(:,:,jk)
378      END DO
379      ! Boundary conditions on the laplacian  (zlu,zlv)
380      CALL lbc_lnk_adj( zlvad, 'V', -1.0_wp )
381      CALL lbc_lnk_adj( zluad, 'U', -1.0_wp )
382      !                                                ! ===============
383      DO jk = 1, jpkm1                                 ! Horizontal slab
384         !                                             ! ===============
385         ! Laplacian
386         ! ---------
387         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
388            DO jj = jpjm1, 2, -1
389               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
390                  zufad   (ji  ,jj  ,jk) = zufad   (ji  ,jj  ,jk) + zlvad(ji  ,jj  ,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
391                  zufad   (ji-1,jj  ,jk) = zufad   (ji-1,jj  ,jk) - zlvad(ji  ,jj  ,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
392                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj)
393                  hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj)
394                  zlvad(ji,jj,jk) = 0.0_wp
395                  zufad   (ji  ,jj  ,jk) = zufad   (ji  ,jj  ,jk) - zluad(ji  ,jj  ,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
396                  zufad   (ji  ,jj-1,jk) = zufad   (ji  ,jj-1,jk) + zluad(ji  ,jj  ,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
397                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zluad(ji,jj,jk) / e1u(ji,jj)
398                  hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) + zluad(ji,jj,jk) / e1u(ji,jj)
399                  zluad(ji,jj,jk) = 0.0_wp
400               END DO
401            END DO
402            rotb_ad(:,:,jk) = rotb_ad(:,:,jk) + zufad(:,:,jk) * fse3f(:,:,jk)
403         ELSE                            ! z-coordinate - full step
404            DO jj = jpjm1, 2, -1
405               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
406                  rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) + zlvad(ji,jj,jk) / e1v(ji,jj)
407                  rotb_ad (ji-1,jj  ,jk) = rotb_ad (ji-1,jj  ,jk) - zlvad(ji,jj,jk) / e1v(ji,jj)
408                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj)
409                  hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj)
410                  zlvad(ji,jj,jk) = 0.0_wp
411                  rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) - zluad(ji,jj,jk) / e2u(ji,jj)
412                  rotb_ad (ji  ,jj-1,jk) = rotb_ad (ji  ,jj-1,jk) + zluad(ji,jj,jk) / e2u(ji,jj)
413                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zluad(ji,jj,jk) / e1u(ji,jj)
414                  hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) + zluad(ji,jj,jk) / e1u(ji,jj)
415                  zlvad(ji,jj,jk) = 0.0_wp
416               END DO
417            END DO
418         ENDIF
419      ENDDO
420      CALL wrk_dealloc( jpi, jpj,      zcuad, zcvad           )
421      CALL wrk_dealloc( jpi, jpj, jpk, zufad, zutad, zluad, zlvad )
422      !
423      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap_adj')
424      !
425   END SUBROUTINE dyn_ldf_bilap_adj
426
427   !!======================================================================
428#endif
429END MODULE dynldf_bilap_tam
Note: See TracBrowser for help on using the repository browser.