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_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynldf_bilap_tam.F90 @ 4575

Last change on this file since 4575 was 4575, checked in by pabouttier, 10 years ago

Cosmetic changes in dynldf_bilap_tam
  • Property svn:executable set to *
File size: 20.0 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      !!   9.0  !  09-12 (F. Vigilant) adjoint of 9.0
269      !!   3.4  !  12-07  (P.-A. Bouttier) 3.4 version
270      !!----------------------------------------------------------------------
271      !! * Arguments
272      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
273      !! * Local declarations
274      INTEGER  ::   ji, jj, jk            ! dummy loop indices
275      REAL(wp) ::   &
276         zuaad, zvaad, zbt, ze2u, ze2v        ! temporary scalars
277      REAL(wp), POINTER, DIMENSION(:,:) ::       &
278         zcuad, zcvad                         ! temporary workspace
279      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
280         zufad, zutad, zluad, zlvad           ! temporary workspace
281      !!----------------------------------------------------------------------
282      !
283      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap_adj')
284      !
285      CALL wrk_alloc( jpi, jpj,      zcuad, zcvad           )
286      CALL wrk_alloc( jpi, jpj, jpk, zufad, zutad, zluad, zlvad )
287      !
288      IF( kt == nitend ) THEN
289         IF(lwp) WRITE(numout,*)
290         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap_adj: bilaplacien operator'
291         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ '
292      ENDIF
293
294      zuaad = 0.0_wp
295      zvaad = 0.0_wp
296
297      zufad(:,:,:) = 0.0_wp
298      zutad(:,:,:) = 0.0_wp
299      zluad(:,:,:) = 0.0_wp
300      zlvad(:,:,:) = 0.0_wp
301
302      zcvad(:,:) = 0.0_wp
303      zcuad(:,:) = 0.0_wp
304
305      DO jk = 1, jpkm1
306
307         ! Bilaplacian
308         ! -----------
309
310         DO jj = jpjm1, 2, -1
311            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
312               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
313               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
314               ! add it to the general momentum trends
315               zvaad = zvaad + va_ad(ji,jj,jk)
316               zuaad = zuaad + ua_ad(ji,jj,jk)
317               ! horizontal biharmonic diffusive trends
318               zufad(ji  ,jj  ,jk) = zufad(ji  ,jj  ,jk) + zvaad / ze2v
319               zufad(ji-1,jj  ,jk) = zufad(ji-1,jj  ,jk) - zvaad / ze2v
320               zutad(ji  ,jj  ,jk) = zutad(ji  ,jj  ,jk) - zvaad / e2v(ji,jj)
321               zutad(ji  ,jj+1,jk) = zutad(ji  ,jj+1,jk) + zvaad / e2v(ji,jj)
322
323               zufad(ji  ,jj  ,jk) = zufad(ji  ,jj  ,jk) - zuaad / ze2u
324               zufad(ji  ,jj-1,jk) = zufad(ji  ,jj-1,jk) + zuaad / ze2u
325               zutad(ji  ,jj  ,jk) = zutad(ji  ,jj  ,jk) - zuaad / e1u(ji,jj)
326               zutad(ji+1,jj  ,jk) = zutad(ji+1,jj  ,jk) + zuaad / e1u(ji,jj)
327               
328               zuaad = 0.0_wp
329               zvaad = 0.0_wp
330            END DO
331         END DO
332
333         !                                             ! ===============
334      END DO                                           !   End of slab
335      !                                                ! ===============
336
337      ! boundary conditions on the laplacian curl and div (zuf,zut)
338!!bug gm no need to do this 2 following lbc...
339      CALL lbc_lnk_adj( zutad, 'T', 1.0_wp )
340      CALL lbc_lnk_adj( zufad, 'F', 1.0_wp )
341
342      DO jk = 1, jpkm1
343
344         ! Third derivative
345         ! ----------------
346
347         ! Laplacian divergence
348         DO jj = jpj, 2, -1
349            DO ji = jpi, fs_2, -1   ! vector opt.
350               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
351               zluad(ji  ,jj  ,jk) = zluad(ji  ,jj  ,jk) + zutad(ji,jj,jk) / zbt
352               zluad(ji-1,jj  ,jk) = zluad(ji-1,jj  ,jk) - zutad(ji,jj,jk) / zbt
353               zlvad(ji  ,jj  ,jk) = zlvad(ji  ,jj  ,jk) + zutad(ji,jj,jk) / zbt
354               zlvad(ji  ,jj-1,jk) = zlvad(ji  ,jj-1,jk) - zutad(ji,jj,jk) / zbt
355
356               zutad(ji,jj,jk) = 0.0_wp
357            END DO
358         END DO
359
360         ! Laplacian Horizontal fluxes
361         DO jj = jpjm1, 1, -1
362            DO ji = fs_jpim1, 1, -1   ! vector opt.
363               zluad(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zluad(ji,jj,jk)
364               zlvad(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlvad(ji,jj,jk)
365            END DO
366         END DO
367
368         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
369         DO jj = jpjm1, 1, -1
370            DO ji = fs_jpim1, 1, -1   ! vector opt.
371               zufad(ji,jj,jk) = fmask(ji,jj,jk) * zufad(ji,jj,jk) * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
372               zcvad(ji  ,jj  ) = zcvad(ji  ,jj  ) - zufad(ji,jj,jk)
373               zcvad(ji+1,jj  ) = zcvad(ji+1,jj  ) + zufad(ji,jj,jk)
374               zcuad(ji  ,jj  ) = zcuad(ji  ,jj  ) + zufad(ji,jj,jk)
375               zcuad(ji  ,jj+1) = zcuad(ji  ,jj+1) - zufad(ji,jj,jk)
376
377               zufad(ji,jj,jk) = 0.0_wp
378            END DO
379         END DO
380
381         ! Contravariant "laplacian"
382         DO jj = 1, jpj
383            DO ji = 1, jpi
384               zlvad(ji,jj,jk) = zlvad(ji,jj,jk) + e2v(ji,jj) * zcvad(ji,jj)
385               zluad(ji,jj,jk) = zluad(ji,jj,jk) + e1u(ji,jj) * zcuad(ji,jj)
386               zcvad(ji,jj) = 0.0_wp
387               zcuad(ji,jj) = 0.0_wp
388            END DO
389         END DO
390
391         ! Multiply by the eddy viscosity coef. (at u- and v-points)
392         zluad(:,:,jk) = zluad(:,:,jk) * fsahmu(:,:,jk)
393         zlvad(:,:,jk) = zlvad(:,:,jk) * fsahmv(:,:,jk)
394
395      END DO
396
397      ! Boundary conditions on the laplacian  (zlu,zlv)
398      CALL lbc_lnk_adj( zlvad, 'V', -1.0_wp )
399      CALL lbc_lnk_adj( zluad, 'U', -1.0_wp )
400
401      !                                                ! ===============
402      DO jk = 1, jpkm1                                 ! Horizontal slab
403         !                                             ! ===============
404         ! Laplacian
405         ! ---------
406
407         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
408            DO jj = jpjm1, 2, -1
409               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
410                  zufad   (ji  ,jj  ,jk) = zufad   (ji  ,jj  ,jk) + zlvad(ji  ,jj  ,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
411                  zufad   (ji-1,jj  ,jk) = zufad   (ji-1,jj  ,jk) - zlvad(ji  ,jj  ,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
412                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj)
413                  hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj)
414                  zlvad(ji,jj,jk) = 0.0_wp
415
416                  zufad   (ji  ,jj  ,jk) = zufad   (ji  ,jj  ,jk) - zluad(ji  ,jj  ,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
417                  zufad   (ji  ,jj-1,jk) = zufad   (ji  ,jj-1,jk) + zluad(ji  ,jj  ,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
418                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zluad(ji,jj,jk) / e1u(ji,jj)
419                  hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) + zluad(ji,jj,jk) / e1u(ji,jj)
420                  zluad(ji,jj,jk) = 0.0_wp
421               END DO
422            END DO
423            rotb_ad(:,:,jk) = rotb_ad(:,:,jk) + zufad(:,:,jk) * fse3f(:,:,jk)
424         ELSE                            ! z-coordinate - full step
425            DO jj = jpjm1, 2, -1
426               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
427                  rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) + zlvad(ji,jj,jk) / e1v(ji,jj)
428                  rotb_ad (ji-1,jj  ,jk) = rotb_ad (ji-1,jj  ,jk) - zlvad(ji,jj,jk) / e1v(ji,jj)
429                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj)
430                  hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj)
431                  zlvad(ji,jj,jk) = 0.0_wp
432                 
433                  rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) - zluad(ji,jj,jk) / e2u(ji,jj)
434                  rotb_ad (ji  ,jj-1,jk) = rotb_ad (ji  ,jj-1,jk) + zluad(ji,jj,jk) / e2u(ji,jj)
435                  hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) - zluad(ji,jj,jk) / e1u(ji,jj)
436                  hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) + zluad(ji,jj,jk) / e1u(ji,jj)
437                  zlvad(ji,jj,jk) = 0.0_wp
438               END DO
439            END DO
440         ENDIF
441      ENDDO
442      CALL wrk_dealloc( jpi, jpj,      zcuad, zcvad           )
443      CALL wrk_dealloc( jpi, jpj, jpk, zufad, zutad, zluad, zlvad )
444      !
445      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap_adj')
446      !
447   END SUBROUTINE dyn_ldf_bilap_adj
448
449   !!======================================================================
450#endif
451END MODULE dynldf_bilap_tam
Note: See TracBrowser for help on using the repository browser.