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.
divcur_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/divcur_tam.F90 @ 4573

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

Fix the case where jpi is odd, in northfold treatment for T-point in div_cur_adj (divcur_tam.F90 module), see Ticket #1270
  • Property svn:executable set to *
File size: 53.2 KB
Line 
1MODULE divcur_tam
2#if defined key_tam
3   !!==============================================================================
4   !!       ***  MODULE  divcur_tam : TANGENT/ADJOINT OF MODULE divcur  ***
5   !!
6   !! Ocean diagnostic variable : horizontal divergence and relative vorticity
7   !!
8   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
9   !!            4.0  ! 1991-11  (G. Madec)
10   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
11   !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
12   !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
13   !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
14   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
15   !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
16   !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
17   !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
18   !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
19   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
20   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here
21   !!==============================================================================
22   !! History of the TAM module:
23   !!             7.0  !  95-01  (F. Van den Berghe)
24   !!             8.0  !  96-04  (A. Weaver)
25   !!             8.1  !  98-02  (A. Weaver)
26   !!             8.2  !  00-08  (A. Weaver)
27   !!             9.0  !  08-06  (A. Vidard) Skeleton
28   !!             9.0  !  08-07  (A. Weaver)
29   !!             9.0  !  08-11  (A. Vidard) Nemo v3
30   !!             9.0  !  09-02  (A. Vidard) cleanup
31   !!             9.0  !  07-12  (P.-A. Bouttier) Nemo v3.4
32   !!----------------------------------------------------------------------
33   !!   div_cur_tan     : Compute the horizontal divergence and relative
34   !!                     vorticity fields (tangent routine)
35   !!   div_cur_adj     : Compute the horizontal divergence and relative
36   !!                     vorticity fields (adjoint routine)
37   !!   div_cur_adj_tst : Test of the adjoint routine
38   !!----------------------------------------------------------------------
39   !! * Modules used
40   USE par_kind
41   USE par_oce
42   USE in_out_manager
43   USE dom_oce
44   USE sbc_oce
45   USE lbclnk
46   USE lbclnk_tam
47   USE oce_tam
48   USE gridrandom
49   USE dotprodfld
50   USE tstool_tam
51   USE lib_mpp         ! MPP library
52   USE wrk_nemo        ! Memory Allocation
53   USE timing          ! Timing
54   USE cla_tam
55   USE sbcrnf_tam
56
57   PRIVATE
58
59   !! * Accessibility
60   PUBLIC div_cur_tan,    &   ! routine called by steptan.F90
61      &   div_cur_adj,    &   ! routine called by stepadj.F90
62      &   div_cur_adj_tst     ! adjoint test routine
63
64   !! * Substitutions
65#  include "domzgr_substitute.h90"
66#  include "vectopt_loop_substitute.h90"
67
68CONTAINS
69
70#if defined key_noslip_accurate
71   !!----------------------------------------------------------------------
72   !!   'key_noslip_accurate'                     2nd order centered scheme
73   !!                                                4th order at the coast
74   !!----------------------------------------------------------------------
75
76   SUBROUTINE div_cur_tan( kt )
77      !!----------------------------------------------------------------------
78      !!                  ***  ROUTINE div_cur_tan  ***
79      !!
80      !! ** Purpose of direct routine :
81      !!      compute the horizontal divergence and the relative
82      !!      vorticity at before and now time-step
83      !!
84      !! ** Method of direct routine :
85      !!      I.  divergence :
86      !!         - save the divergence computed at the previous time-step
87      !!      (note that the Asselin filter has not been applied on hdivb)
88      !!         - compute the now divergence given by :
89      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
90      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
91      !!      above expression
92      !!         - apply lateral boundary conditions on hdivn
93      !!      II. vorticity :
94      !!         - save the curl computed at the previous time-step
95      !!            rotb = rotn
96      !!      (note that the Asselin time filter has not been applied to rotb)
97      !!         - compute the now curl in tensorial formalism:
98      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
99      !!         - apply lateral boundary conditions on rotn through a call
100      !!      of lbc_lnk routine.
101      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
102      !!      the no-slip boundary condition is computed using Schchepetkin
103      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
104      !!      For example, along east coast, the one-sided finite difference
105      !!      approximation used for di[v] is:
106      !!         di[e2v vn] =  1/(e1f*e2f)
107      !!                    * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
108      !!
109      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
110      !!              - update rotb , rotn , the before & now rel. vorticity
111      !!
112      !! History of the direct routine:
113      !!   8.2  !  00-03  (G. Madec)  no slip accurate
114      !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
115      !!        !  05-01  (J. Chanut, A. Sellar) unstructured open boundaries
116      !! History of the TAM routine:
117      !!   9.0  !  08-06  (A. Vidard) Skeleton
118      !!        !  08-07  (A. Weaver)
119      !!        !  08-11  (A. Vidard) Nemo v3
120      !!----------------------------------------------------------------------
121      !! * Arguments
122      INTEGER, INTENT( in ) :: &
123         & kt     ! ocean time-step index
124
125      !! * Local declarations
126      INTEGER :: &
127         & ji,  & ! dummy loop indices
128         & jj,  &
129         & jk
130      INTEGER :: &
131         & ii,  & ! temporary integer
132         & ij,  &
133         & jl,  &
134         & ijt, &
135         & iju
136      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
137         & zwu    ! Workspace
138      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
139         & zwv    ! Workspace
140      !!----------------------------------------------------------------------
141      IF( nn_timing == 1 )  CALL timing_start('div_cur_tan')
142      !
143      CALL wrk_alloc( jpi  , jpj+2, zwu               )
144      CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 )
145      !
146      IF( kt == nit000 ) THEN
147         IF(lwp) WRITE(numout,*)
148         IF(lwp) WRITE(numout,*) 'div_cur_tan : horizontal velocity', &
149            &                    ' divergence and relative vorticity'
150         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   NOT optimal for auto-tasking case'
151      ENDIF
152
153      !                                                ! ===============
154      DO jk = 1, jpkm1                                 ! Horizontal slab
155         !                                             ! ===============
156
157         hdivb_tl(:,:,jk) = hdivn_tl(:,:,jk)    ! time swap of div arrays
158         rotb_tl (:,:,jk) = rotn_tl (:,:,jk)    ! time swap of rot arrays
159
160         !                                             ! --------
161         ! Horizontal divergence                       !   div
162         !                                             ! --------
163         DO jj = 2, jpjm1
164            DO ji = fs_2, fs_jpim1   ! vector opt.
165               hdivn_tl(ji,jj,jk) =   &
166                  & (   e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un_tl(ji  ,jj  ,jk) &
167                  &   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un_tl(ji-1,jj  ,jk) &
168                  &   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn_tl(ji  ,jj  ,jk) &
169                  &   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn_tl(ji  ,jj-1,jk) &
170                  & ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
171            END DO
172         END DO
173         !
174         IF ( .NOT. Agrif_Root() ) then
175            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn_tl(nlci-1 , :     ,jk) = 0.e0      ! east
176            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn_tl(2      , :     ,jk) = 0.e0      ! west
177            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn_tl(:      ,nlcj-1 ,jk) = 0.e0      ! north
178            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn_tl(:      ,2      ,jk) = 0.e0      ! south
179         ENDIF
180         !
181         !                                             ! --------
182         ! relative vorticity                          !   rot
183         !                                             ! --------
184         ! contravariant velocity (extended for lateral b.c.)
185         ! inside the model domain
186         DO jj = 1, jpj
187            DO ji = 1, jpi
188               zwu(ji,jj) = e1u(ji,jj) * un_tl(ji,jj,jk)
189               zwv(ji,jj) = e2v(ji,jj) * vn_tl(ji,jj,jk)
190            END DO
191         END DO
192
193         ! East-West boundary conditions
194         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
195            zwv(  0  ,:) = zwv(jpi-2,:)
196            zwv( -1  ,:) = zwv(jpi-3,:)
197            zwv(jpi+1,:) = zwv(  3  ,:)
198            zwv(jpi+2,:) = zwv(  4  ,:)
199         ELSE
200            zwv(  0  ,:) = 0.0_wp
201            zwv( -1  ,:) = 0.0_wp
202            zwv(jpi+1,:) = 0.0_wp
203            zwv(jpi+2,:) = 0.0_wp
204         ENDIF
205
206         ! North-South boundary conditions
207         IF( nperio == 3 .OR. nperio == 4 ) THEN
208            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre
209            zwu(jpi,jpj+1) = 0.0_wp
210            zwu(jpi,jpj+2) = 0.0_wp
211            DO ji = 1, jpi-1
212               iju = jpi - ji + 1
213               zwu(ji,jpj+1) = - zwu(iju,jpj-3)
214               zwu(ji,jpj+2) = - zwu(iju,jpj-4)
215            END DO
216         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
217            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre
218            zwu(jpi,jpj+1) = 0.0_wp
219            zwu(jpi,jpj+2) = 0.0_wp
220            DO ji = 1, jpi-1
221               iju = jpi - ji
222               zwu(ji,jpj  ) = - zwu(iju,jpj-1)
223               zwu(ji,jpj+1) = - zwu(iju,jpj-2)
224               zwu(ji,jpj+2) = - zwu(iju,jpj-3)
225            END DO
226            DO ji = -1, jpi+2
227               ijt = jpi - ji + 1
228               zwv(ji,jpj) = - zwv(ijt,jpj-2)
229            END DO
230            DO ji = jpi/2+1, jpi+2
231               ijt = jpi - ji + 1
232               zwv(ji,jpjm1) = - zwv(ijt,jpjm1)
233            END DO
234         ELSE
235            ! closed
236            zwu(:,jpj+1) = 0.0_wp
237            zwu(:,jpj+2) = 0.0_wp
238         ENDIF
239
240         ! relative vorticity (vertical component of the velocity curl)
241         DO jj = 1, jpjm1
242            DO ji = 1, fs_jpim1   ! vector opt.
243               rotn_tl(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
244                  &                 - zwu(ji  ,jj+1) + zwu(ji,jj)  )   &
245                  &         * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
246            END DO
247         END DO
248
249         ! second order accurate scheme along straight coast
250         DO jl = 1, npcoa(1,jk)
251            ii = nicoa(jl,1,jk)
252            ij = njcoa(jl,1,jk)
253            rotn_tl(ii,ij,jk) =  1.0_wp / ( e1f(ii,ij) * e2f(ii,ij) )  &
254               & * ( + 4.0_wp * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2_wp * zwv(ii+3,ij) )
255         END DO
256         DO jl = 1, npcoa(2,jk)
257            ii = nicoa(jl,2,jk)
258            ij = njcoa(jl,2,jk)
259            rotn_tl(ii,ij,jk) =  1.0_wp / ( e1f(ii,ij) * e2f(ii,ij) )  &
260               & * ( - 4.0_wp * zwv(ii,ij)   + zwv(ii-1,ij) - 0.2_wp * zwv(ii-2,ij) )
261         END DO
262         DO jl = 1, npcoa(3,jk)
263            ii = nicoa(jl,3,jk)
264            ij = njcoa(jl,3,jk)
265            rotn_tl(ii,ij,jk) = -1.0_wp / ( e1f(ii,ij) * e2f(ii,ij) )  &
266               & * ( + 4.0_wp * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2_wp * zwu(ii,ij+3) )
267         END DO
268         DO jl = 1, npcoa(4,jk)
269            ii = nicoa(jl,4,jk)
270            ij = njcoa(jl,4,jk)
271            rotn_tl(ii,ij,jk) = -1.0_wp / ( e1f(ii,ij) * e2f(ii,ij) )  &
272               & * ( -4.0_wp * zwu(ii,ij)    + zwu(ii,ij-1) - 0.2_wp * zwu(ii,ij-2) )
273         END DO
274
275         !                                             ! ===============
276      END DO                                           !   End of slab
277      !                                                ! ===============
278      IF( ln_rnf      )   CALL sbc_rnf_div_tan( hdivn_tl )       ! runoffs (update hdivn field)
279      IF( nn_cla == 1 )   CALL cla_div_tan    ( kt )             ! Cross Land Advection (Update Hor. divergence)
280
281      ! 4. Lateral boundary conditions on hdivn and rotn
282      ! ---------------------------------=======---======
283      CALL lbc_lnk( hdivn_tl, 'T', 1.0_wp )   ! T-point, no sign change
284      CALL lbc_lnk( rotn_tl , 'F', 1.0_wp )   ! F-point, no sign change
285      !
286      CALL wrk_dealloc( jpi  , jpj+2, zwu               )
287      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 )
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('div_cur_tan')
290   END SUBROUTINE div_cur_tan
291
292   SUBROUTINE div_cur_adj( kt )
293      !!----------------------------------------------------------------------
294      !!                  ***  ROUTINE div_cur_adj  ***
295      !!
296      !! ** Purpose of direct routine :
297      !!      compute the horizontal divergence and the relative
298      !!      vorticity at before and now time-step
299      !!
300      !! ** Method of direct routine :
301      !!      I.  divergence :
302      !!         - save the divergence computed at the previous time-step
303      !!      (note that the Asselin filter has not been applied on hdivb)
304      !!         - compute the now divergence given by :
305      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
306      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
307      !!      above expression
308      !!         - apply lateral boundary conditions on hdivn
309      !!      II. vorticity :
310      !!         - save the curl computed at the previous time-step
311      !!            rotb = rotn
312      !!      (note that the Asselin time filter has not been applied to rotb)
313      !!         - compute the now curl in tensorial formalism:
314      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
315      !!         - apply lateral boundary conditions on rotn through a call
316      !!      of lbc_lnk routine.
317      !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
318      !!      the no-slip boundary condition is computed using Schchepetkin
319      !!      and O'Brien (1996) scheme (i.e. 4th order at the coast).
320      !!      For example, along east coast, the one-sided finite difference
321      !!      approximation used for di[v] is:
322      !!         di[e2v vn] =  1/(e1f*e2f)
323      !!                    * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
324      !!
325      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
326      !!              - update rotb , rotn , the before & now rel. vorticity
327      !!
328      !! History of the direct routine:
329      !!   8.2  !  00-03  (G. Madec)  no slip accurate
330      !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
331      !! History of the TAM routine:
332      !!   9.0  !  08-06  (A. Vidard) Skeleton
333      !!   9.0  !  08-07  (A. Weaver)
334      !!----------------------------------------------------------------------
335      !! * Arguments
336      INTEGER, INTENT( in ) :: &
337         & kt     ! ocean time-step index
338
339      !! * Local declarations
340      INTEGER :: &
341         & ji,  & ! dummy loop indices
342         & jj,  &
343         & jk
344      INTEGER :: &
345         & ii,  & ! temporary integer
346         & ij,  &
347         & jl,  &
348         & ijt, &
349         & iju
350      REAL(wp) :: &
351         & zdiv, &
352         & zdju
353      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
354         & zwu    ! Workspace
355      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
356         & zwv    ! Workspace
357      REAL(wp) :: &
358         & ztmp
359      !!----------------------------------------------------------------------
360      IF( nn_timing == 1 )  CALL timing_start('div_cur_adj')
361      !
362      CALL wrk_alloc( jpi  , jpj+2, zwu               )
363      CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 )
364      !
365      IF( kt == nitend ) THEN
366         IF(lwp) WRITE(numout,*)
367         IF(lwp) WRITE(numout,*) 'div_cur_adj : horizontal velocity', &
368            &                    ' divergence and relative vorticity'
369         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   NOT optimal for auto-tasking case'
370      ENDIF
371
372      ! 4. Lateral boundary conditions on hdivn and rotn
373      ! ---------------------------------=======---======
374      CALL lbc_lnk_adj( rotn_ad , 'F', 1.0_wp )   ! F-point, no sign change
375      CALL lbc_lnk_adj( hdivn_ad, 'T', 1.0_wp )   ! T-point, no sign change
376
377      IF( nn_cla == 1 )   CALL cla_div_adj    ( kt )             ! Cross Land Advection (Update Hor. divergence)
378      IF( ln_rnf      )   CALL sbc_rnf_div_adj( hdivn_ad )       ! runoffs (update hdivn field)
379      !                                                ! ===============
380      DO jk = jpkm1, 1, -1                             ! Horizontal slab
381         !                                             ! ===============
382         ! local adjoint workspace initialization
383         zwu(:,:) = 0.0_wp
384         zwv(:,:) = 0.0_wp
385         !                                             ! --------
386         ! relative vorticity                          !   rot
387         !                                             ! --------
388         DO jl = npcoa(4,jk), 1, -1
389            ii = nicoa(jl,4,jk)
390            ij = njcoa(jl,4,jk)
391            rotn_ad(ii,ij,jk) = -1.0_wp * rotn_ad(ii,ij,jk) &
392               &                / ( e1f(ii,ij) * e2f(ii,ij) )
393            zwu(ii,ij  ) = zwu(ii,ij  ) - 4.0_wp * rotn_ad(ii,ij,jk)
394            zwu(ii,ij-1) = zwu(ii,ij-1) +          rotn_ad(ii,ij,jk)
395            zwu(ii,ij-2) = zwu(ii,ij-2) - 0.2_wp * rotn_ad(ii,ij,jk)
396            rotn_ad(ii,ij,jk) = 0.0_wp
397         END DO
398         DO jl = npcoa(3,jk), 1, -1
399            ii = nicoa(jl,3,jk)
400            ij = njcoa(jl,3,jk)
401            rotn_ad(ii,ij,jk) = -1.0_wp * rotn_ad(ii,ij,jk) &
402               &                / ( e1f(ii,ij) * e2f(ii,ij) )
403            zwu(ii,ij+1) = zwu(ii,ij+1) + 4.0_wp * rotn_ad(ii,ij,jk)
404            zwu(ii,ij+2) = zwu(ii,ij+2) -          rotn_ad(ii,ij,jk)
405            zwu(ii,ij+3) = zwu(ii,ij+3) + 0.2_wp * rotn_ad(ii,ij,jk)
406            rotn_ad(ii,ij,jk) = 0.0_wp
407         END DO
408         DO jl = npcoa(2,jk), 1, -1
409            ii = nicoa(jl,2,jk)
410            ij = njcoa(jl,2,jk)
411            rotn_ad(ii,ij,jk) =  1.0_wp * rotn_ad(ii,ij,jk) &
412               &                / ( e1f(ii,ij) * e2f(ii,ij) )
413            zwv(ii  ,ij) = zwv(ii  ,ij) - 4.0_wp * rotn_ad(ii,ij,jk)
414            zwv(ii-1,ij) = zwv(ii-1,ij) +          rotn_ad(ii,ij,jk)
415            zwv(ii-2,ij) = zwv(ii-2,ij) - 0.2_wp * rotn_ad(ii,ij,jk)
416            rotn_ad(ii,ij,jk) = 0.0_wp
417         END DO
418         ! second order accurate scheme along straight coast
419         DO jl = npcoa(1,jk), 1, -1
420            ii = nicoa(jl,1,jk)
421            ij = njcoa(jl,1,jk)
422            rotn_ad(ii,ij,jk) =  1.0_wp * rotn_ad(ii,ij,jk) &
423               &                / ( e1f(ii,ij) * e2f(ii,ij) )
424            zwv(ii+1,ij) = zwv(ii+1,ij) + 4.0_wp * rotn_ad(ii,ij,jk)
425            zwv(ii+2,ij) = zwv(ii+2,ij) -          rotn_ad(ii,ij,jk)
426            zwv(ii+3,ij) = zwv(ii+3,ij) + 0.2_wp * rotn_ad(ii,ij,jk)
427            rotn_ad(ii,ij,jk) = 0.0_wp
428         END DO
429         ! relative vorticity (vertical component of the velocity curl)
430         DO jj = jpjm1, 1, -1
431            DO ji = fs_jpim1, 1, -1   ! vector opt.
432               rotn_ad(ji,jj,jk) = rotn_ad(ji,jj,jk) * fmask(ji,jj,jk) &
433                  &                / ( e1f(ji,jj) * e2f(ji,jj) )
434               zwv(ji  ,jj  ) = zwv(ji  ,jj  ) - rotn_ad(ji,jj,jk)
435               zwu(ji  ,jj  ) = zwu(ji  ,jj  ) + rotn_ad(ji,jj,jk)
436               zwu(ji  ,jj+1) = zwu(ji  ,jj+1) - rotn_ad(ji,jj,jk)
437               zwv(ji+1,jj  ) = zwv(ji+1,jj  ) + rotn_ad(ji,jj,jk)
438               rotn_ad(ji,jj,jk) = 0.0
439            END DO
440         END DO
441         ! North-South boundary conditions
442         IF( nperio == 3 .OR. nperio == 4 ) THEN
443            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degree
444            DO ji = jpi-1, 1, -1
445               iju = jpi - ji + 1
446               zwu(iju,jpj-4) = zwu(iju,jpj-4) - zwu(ji,jpj+2)
447               zwu(ji ,jpj+2) = 0.0_wp
448               zwu(iju,jpj-3) = zwu(iju,jpj-3) - zwu(ji,jpj+1)
449               zwu(ji ,jpj+1) = 0.0_wp
450            END DO
451            zwu(jpi,jpj+2) = 0.0_wp
452            zwu(jpi,jpj+1) = 0.0_wp
453         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
454            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degree
455            DO ji = jpi+2, jpi/2+1, -1 ! if jpi is odd, you can get ji=ijt hence the use of ztmp
456               ijt = jpi - ji + 1
457               ztmp = zwv(ji, jpjm1)
458               zwv(ji, jpjm1) = 0.0_wp
459               zwv(ijt,jpjm1) = zwv(ijt,jpjm1) - ztmp
460            END DO
461            DO ji = jpi+2, -1, -1
462               ijt = jpi - ji + 1
463               zwv(ijt,jpj-2) = zwv(ijt,jpj-2) - zwv(ji,jpj  )
464               zwv(ji ,jpj  ) = 0.0_wp
465            END DO
466            DO ji = jpi-1, 1, -1
467               iju = jpi - ji
468               zwu(iju,jpj-3) = zwu(iju,jpj-3) - zwu(ji,jpj+2)
469               zwu(ji ,jpj+2) = 0.0_wp
470               zwu(iju,jpj-2) = zwu(iju,jpj-2) - zwu(ji,jpj+1)
471               zwu(ji ,jpj+1) = 0.0_wp
472               zwu(iju,jpj-1) = zwu(iju,jpj-1) - zwu(ji,jpj  )
473               zwu(ji ,jpj  ) = 0.0_wp
474            END DO
475            zwu(jpi,jpj+2) = 0.0_wp
476            zwu(jpi,jpj+1) = 0.0_wp
477         ELSE
478            ! closed
479            zwu(:,jpj+2) = 0.0_wp
480            zwu(:,jpj+1) = 0.0_wp
481         ENDIF
482         ! East-West boundary conditions
483         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
484            zwv(  4  ,:) = zwv(  4  ,:) + zwv(jpi+2,:)
485            zwv(jpi+2,:) = 0.0_wp
486            zwv(  3  ,:) = zwv(  3  ,:) + zwv(jpi+1,:)
487            zwv(jpi+1,:) = 0.0_wp
488            zwv(jpi-3,:) = zwv(jpi-3,:) + zwv( -1  ,:)
489            zwv( -1  ,:) = 0.0_wp
490            zwv(jpi-2,:) = zwv(jpi-2,:) + zwv(  0  ,:)
491            zwv(  0  ,:) = 0.0_wp
492         ELSE
493            zwv(jpi+2,:) = 0.0_wp
494            zwv(jpi+1,:) = 0.0_wp
495            zwv( -1  ,:) = 0.0_wp
496            zwv(  0  ,:) = 0.0_wp
497         ENDIF
498         ! contravariant velocity (extended for lateral b.c.)
499         ! inside the model domain
500         DO jj = jpj, 1, -1
501            DO ji = jpi, 1, -1
502               vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e2v(ji,jj) * zwv(ji,jj)
503               un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e1u(ji,jj) * zwu(ji,jj)
504            END DO
505         END DO
506         !
507         IF( .NOT. AGRIF_Root() ) THEN
508            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn_ad(nlci-1 , :     ,jk) = 0.0_wp      ! east
509            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn_ad(2      , :     ,jk) = 0.0_wp      ! west
510            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn_ad(:      ,nlcj-1 ,jk) = 0.0_wp      ! north
511            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn_ad(:      ,2      ,jk) = 0.0_wp      ! south
512         ENDIF
513         !                                             ! --------
514         ! Horizontal divergence                       !   div
515         !                                             ! --------
516         DO jj = jpjm1, 2, -1
517            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
518               hdivn_ad(ji,jj,jk) =  hdivn_ad(ji,jj,jk)  &
519                  &                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
520               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
521                  &                  + e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) &
522                  &                                   * hdivn_ad(ji,jj,jk)
523               un_ad(ji-1,jj  ,jk) = un_ad(ji-1,jj  ,jk) &
524                  &                  - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) &
525                  &                                   * hdivn_ad(ji,jj,jk)
526               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
527                  &                  + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) &
528                  &                                   * hdivn_ad(ji,jj,jk)
529               vn_ad(ji  ,jj-1,jk) = vn_ad(ji  ,jj-1,jk) &
530                  &                  - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) &
531                  &                                   * hdivn_ad(ji,jj,jk)
532               hdivn_ad(ji,jj,jk) = 0.0_wp
533            END DO
534         END DO
535         rotn_ad (:,:,jk) = rotn_ad (:,:,jk) + rotb_ad (:,:,jk)  ! time swap
536         rotb_ad (:,:,jk) = 0.0_wp
537         hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + hdivb_ad(:,:,jk)  ! time swap
538         hdivb_ad(:,:,jk) = 0.0_wp
539         !                                             ! ===============
540      END DO                                           !   End of slab
541      !                                                ! ===============
542      CALL wrk_dealloc( jpi  , jpj+2, zwu               )
543      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 )
544      !
545      IF( nn_timing == 1 )  CALL timing_stop('div_cur_adj')
546
547   END SUBROUTINE div_cur_adj
548
549#else
550   !!----------------------------------------------------------------------
551   !!   Default option                           2nd order centered schemes
552   !!----------------------------------------------------------------------
553   SUBROUTINE div_cur_tan( kt )
554      !!----------------------------------------------------------------------
555      !!                  ***  ROUTINE div_cur_tan  ***
556      !!
557      !! ** Purpose of direct routine :
558      !!      compute the horizontal divergence and the relative
559      !!      vorticity at before and now time-step
560      !!
561      !! ** Method of direct routine :
562      !!              - Divergence:
563      !!      - save the divergence computed at the previous time-step
564      !!      (note that the Asselin filter has not been applied on hdivb)
565      !!      - compute the now divergence given by :
566      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
567      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
568      !!      above expression
569      !!      - apply lateral boundary conditions on hdivn
570      !!              - Relavtive Vorticity :
571      !!      - save the curl computed at the previous time-step (rotb = rotn)
572      !!      (note that the Asselin time filter has not been applied to rotb)
573      !!      - compute the now curl in tensorial formalism:
574      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
575      !!      - apply lateral boundary conditions on rotn through a call to
576      !!      routine lbc_lnk routine.
577      !!      Note: Coastal boundary condition: lateral friction set through
578      !!      the value of fmask along the coast (see dommsk.F90) and shlat
579      !!      (namelist parameter)
580      !!
581      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
582      !!              - update rotb , rotn , the before & now rel. vorticity
583      !!
584      !! History of the direct routine:
585      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
586      !!   4.0  !  91-11  (G. Madec)
587      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
588      !!   7.0  !  96-01  (G. Madec)  s-coordinates
589      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
590      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
591      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
592      !!        !  05-01  (J. Chanut) Unstructured open boundaries
593      !! History of the TAM routine:
594      !!   9.0  !  08-06  (A. Vidard) Skeleton
595      !!        !  08-07  (A. Weaver) tangent of the 02-09 version
596      !!        !  08-11  (A. Vidard) tangent of the 05-01 version
597      !!----------------------------------------------------------------------
598      !! * Arguments
599      INTEGER, INTENT( in ) :: &
600         & kt     ! ocean time-step index
601
602      !! * Local declarations
603      INTEGER :: &
604         & ji,  & ! dummy loop indices
605         & jj,  &
606         & jk
607      !!----------------------------------------------------------------------
608      IF( nn_timing == 1 )  CALL timing_start('div_cur_tan')
609      !
610      IF( kt == nit000 ) THEN
611         IF(lwp) WRITE(numout,*)
612         IF(lwp) WRITE(numout,*) 'div_cur_tan : horizontal velocity divergence and'
613         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   relative vorticity'
614      ENDIF
615      !                                                ! ===============
616      DO jk = 1, jpkm1                                 ! Horizontal slab
617         !                                             ! ===============
618         hdivb_tl(:,:,jk) = hdivn_tl(:,:,jk)    ! time swap of div arrays
619         rotb_tl (:,:,jk) = rotn_tl (:,:,jk)    ! time swap of rot arrays
620         !                                             ! --------
621         ! Horizontal divergence                       !   div
622         !                                             ! --------
623         DO jj = 2, jpjm1
624            DO ji = fs_2, fs_jpim1   ! vector opt.
625               hdivn_tl(ji,jj,jk) =   &
626                  & (   e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un_tl(ji  ,jj  ,jk) &
627                  &   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un_tl(ji-1,jj  ,jk) &
628                  &   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn_tl(ji  ,jj  ,jk) &
629                  &   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn_tl(ji  ,jj-1,jk) &
630                  & ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
631            END DO
632         END DO
633         !                                             ! --------
634         ! relative vorticity                          !   rot
635         !                                             ! --------
636         DO jj = 1, jpjm1
637            DO ji = 1, fs_jpim1   ! vector opt.
638               rotn_tl(ji,jj,jk) = (   e2v(ji+1,jj  ) * vn_tl(ji+1,jj  ,jk) &
639                  &                  - e2v(ji  ,jj  ) * vn_tl(ji  ,jj  ,jk) &
640                  &                  - e1u(ji  ,jj+1) * un_tl(ji  ,jj+1,jk) &
641                  &                  + e1u(ji  ,jj  ) * un_tl(ji  ,jj  ,jk) &
642                  &                ) * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
643            END DO
644         END DO
645         !                                             ! ===============
646      END DO                                           !   End of slab
647      !                                                ! ===============
648      IF( ln_rnf      )   CALL sbc_rnf_div_tan( hdivn_tl )       ! runoffs (update hdivn field)
649      IF( nn_cla == 1 )   CALL cla_div_tan    ( kt )             ! Cross Land Advection (update hdivn field)
650      !!
651      CALL lbc_lnk( hdivn_tl, 'T', 1. )
652      CALL lbc_lnk( rotn_tl , 'F', 1. )     ! lateral boundary cond. (no sign change)
653      !
654      IF( nn_timing == 1 )  CALL timing_stop('div_cur_tan')
655   END SUBROUTINE div_cur_tan
656
657   SUBROUTINE div_cur_adj( kt )
658      !!----------------------------------------------------------------------
659      !!                  ***  ROUTINE div_cur_adj  ***
660      !!
661      !! ** Purpose of direct routine :
662      !!      compute the horizontal divergence and the relative
663      !!      vorticity at before and now time-step
664      !!
665      !! ** Method of direct routine :
666      !!              - Divergence:
667      !!      - save the divergence computed at the previous time-step
668      !!      (note that the Asselin filter has not been applied on hdivb)
669      !!      - compute the now divergence given by :
670      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
671      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
672      !!      above expression
673      !!      - apply lateral boundary conditions on hdivn
674      !!              - Relavtive Vorticity :
675      !!      - save the curl computed at the previous time-step (rotb = rotn)
676      !!      (note that the Asselin time filter has not been applied to rotb)
677      !!      - compute the now curl in tensorial formalism:
678      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
679      !!      - apply lateral boundary conditions on rotn through a call to
680      !!      routine lbc_lnk routine.
681      !!      Note: Coastal boundary condition: lateral friction set through
682      !!      the value of fmask along the coast (see dommsk.F90) and shlat
683      !!      (namelist parameter)
684      !!
685      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
686      !!              - update rotb , rotn , the before & now rel. vorticity
687      !!
688      !! History of the direct routine:
689      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
690      !!   4.0  !  91-11  (G. Madec)
691      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
692      !!   7.0  !  96-01  (G. Madec)  s-coordinates
693      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
694      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
695      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
696      !! History of the TAM routine:
697      !!   9.0  !  08-06  (A. Vidard) Skeleton
698      !!   9.0  !  08-07  (A. Weaver)
699      !!----------------------------------------------------------------------
700      !! * Arguments
701      INTEGER, INTENT( in ) :: &
702         & kt     ! ocean time-step index
703
704      !! * Local declarations
705      INTEGER :: &
706         & ji,  & ! dummy loop indices
707         & jj,  &
708         & jk
709      !!----------------------------------------------------------------------
710      !
711      if( nn_timing == 1 )  call timing_start('div_cur_adj')
712      !
713      IF( kt == nitend ) THEN
714         IF(lwp) WRITE(numout,*)
715         IF(lwp) WRITE(numout,*) 'div_cur_adj : horizontal velocity divergence and'
716         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   relative vorticity'
717      ENDIF
718      ! 4. Lateral boundary conditions on hdivn and rotn
719      ! ---------------------------------=======---======
720      CALL lbc_lnk_adj( rotn_ad , 'F', 1.0_wp )     ! F-point, no sign change
721      CALL lbc_lnk_adj( hdivn_ad, 'T', 1.0_wp )     ! T-point, no sign change
722      !!
723      IF( nn_cla == 1 )   CALL cla_div_adj    ( kt )             ! Cross Land Advection (update hdivn field)
724      IF( ln_rnf      )   CALL sbc_rnf_div_adj( hdivn_ad )       ! runoffs (update hdivn field)
725      !                                                ! ===============
726      DO jk = jpkm1, 1, -1                             ! Horizontal slab
727         !                                             ! ===============
728         !                                             ! --------
729         ! relative vorticity                          !   rot
730         !                                             ! --------
731         DO jj = jpjm1, 1, -1
732            DO ji = fs_jpim1, 1, -1   ! vector opt.
733               rotn_ad(ji,jj,jk) =  rotn_ad(ji,jj,jk) * fmask(ji,jj,jk) &
734                  &                       / ( e1f(ji,jj) * e2f(ji,jj) )
735               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
736                  &                  + e1u(ji  ,jj  ) * rotn_ad(ji,jj,jk)
737               un_ad(ji  ,jj+1,jk) = un_ad(ji  ,jj+1,jk) &
738                  &                  - e1u(ji  ,jj+1) * rotn_ad(ji,jj,jk)
739               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
740                  &                  - e2v(ji  ,jj  ) * rotn_ad(ji,jj,jk)
741               vn_ad(ji+1,jj  ,jk) = vn_ad(ji+1,jj  ,jk) &
742                  &                  + e2v(ji+1,jj  ) * rotn_ad(ji,jj,jk)
743               rotn_ad(ji,jj,jk) = 0.0_wp
744            END DO
745         END DO
746         !                                             ! --------
747         ! Horizontal divergence                       !   div
748         !                                             ! --------
749         DO jj = jpjm1, 2, -1
750            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
751               hdivn_ad(ji,jj,jk) =  hdivn_ad(ji,jj,jk)  &
752                  &                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
753               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
754                  &                  + e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) &
755                  &                                   * hdivn_ad(ji,jj,jk)
756               un_ad(ji-1,jj  ,jk) = un_ad(ji-1,jj  ,jk) &
757                  &                  - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) &
758                  &                                   * hdivn_ad(ji,jj,jk)
759               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
760                  &                  + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) &
761                  &                                   * hdivn_ad(ji,jj,jk)
762               vn_ad(ji  ,jj-1,jk) = vn_ad(ji  ,jj-1,jk) &
763                  &                  - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) &
764                  &                                   * hdivn_ad(ji,jj,jk)
765               hdivn_ad(ji,jj,jk) = 0.0_wp
766            END DO
767         END DO
768         !
769         rotn_ad (:,:,jk) = rotn_ad (:,:,jk) + rotb_ad (:,:,jk)  ! time swap
770         rotb_ad (:,:,jk) = 0.0_wp
771         hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + hdivb_ad(:,:,jk)  ! time swap
772         hdivb_ad(:,:,jk) = 0.0_wp
773         !                                             ! ===============
774      END DO                                           !   End of slab
775      !                                                ! ===============
776      if( nn_timing == 1 )  call timing_stop('div_cur_adj')
777      !
778   END SUBROUTINE div_cur_adj
779
780#endif
781
782   SUBROUTINE div_cur_adj_tst( kumadt )
783      !!-----------------------------------------------------------------------
784      !!
785      !!          ***  ROUTINE div_cur_adj_tst : TEST OF div_cur_adj  ***
786      !!
787      !! ** Purpose : Test the adjoint routine.
788      !!
789      !! ** Method  : Verify the scalar product
790      !!
791      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
792      !!
793      !!              where  L   = tangent routine
794      !!                     L^T = adjoint routine
795      !!                     W   = diagonal matrix of scale factors
796      !!                     dx  = input perturbation (random field)
797      !!                     dy  = L dx
798      !!
799      !! ** Action  : Separate tests are applied for the following dx and dy:
800      !!
801      !!      1) dx = ( un_tl, vn_tl ) and
802      !!         dy = ( hdivn_tl )
803      !!      2) dx = ( un_tl, vn_tl ) and
804      !!         dy = ( rotntl )
805      !!
806      !! History :
807      !!        ! 08-07 (A. Weaver)
808      !!-----------------------------------------------------------------------
809
810      !! * Modules used
811      !! * Arguments
812      INTEGER, INTENT(IN) :: &
813         & kumadt             ! Output unit
814
815      INTEGER :: &
816         & ji,    &        ! dummy loop indices
817         & jj,    &
818         & jk
819
820      !! * Local declarations
821      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
822         & zun_tlin,     & ! Tangent input: now u-velocity
823         & zvn_tlin,     & ! Tangent input: now v-velocity
824         & zhdivn_tlin,  & ! Tangent input: now horizontal divergence
825         & zrotn_tlin,   & ! Tangent input: now relative vorticity
826         & zhdivb_tlout, & ! Tangent output: before horizontal divergence
827         & zhdivn_tlout, & ! Tangent output: now horizontal divergence
828         & zrotb_tlout,  & ! Tangent output: before relative vorticity
829         & zrotn_tlout,  & ! Tangent output: now relative vorticity
830         & zhdivb_adin,  & ! Adjoint input: before horizontal divergence
831         & zhdivn_adin,  & ! Adjoint input: now horizontal divergence
832         & zrotb_adin,   & ! Adjoint input: before relative vorticity
833         & zrotn_adin,   & ! Adjoint input: now relative vorticity
834         & zun_adout,    & ! Adjoint output: now u-velocity
835         & zvn_adout,    & ! Adjoint output: now v-velocity
836         & zhdivn_adout, & ! Adjoint output: now horizontal divergence
837         & zrotn_adout,  & ! Adjoint output: now relative vorticity
838         & znu,          & ! 3D random field for u
839         & znv             ! 3D random field for v
840
841      REAL(KIND=wp) :: &
842                           ! random field standard deviation for:
843         & zsp1,         & ! scalar product involving the tangent routine
844         & zsp1_1,       & !   scalar product components
845         & zsp1_2,       &
846         & zsp1_3,       & !
847         & zsp1_4,       &
848         & zsp2,         & ! scalar product involving the adjoint routine
849         & zsp2_1,       & !   scalar product components
850         & zsp2_2,       &
851         & zsp2_3,       &
852         & zsp2_4
853
854      CHARACTER(LEN=14) :: cl_name
855
856      ! Allocate memory
857
858      ALLOCATE( &
859         & zun_tlin(jpi,jpj,jpk),     &
860         & zvn_tlin(jpi,jpj,jpk),     &
861         & zhdivn_tlin(jpi,jpj,jpk),  &
862         & zrotn_tlin(jpi,jpj,jpk),   &
863         & zhdivb_tlout(jpi,jpj,jpk), &
864         & zhdivn_tlout(jpi,jpj,jpk), &
865         & zrotb_tlout(jpi,jpj,jpk),  &
866         & zrotn_tlout(jpi,jpj,jpk),  &
867         & zhdivb_adin(jpi,jpj,jpk),  &
868         & zhdivn_adin(jpi,jpj,jpk),  &
869         & zrotb_adin(jpi,jpj,jpk),   &
870         & zrotn_adin(jpi,jpj,jpk),   &
871         & zun_adout(jpi,jpj,jpk),    &
872         & zvn_adout(jpi,jpj,jpk),    &
873         & zhdivn_adout(jpi,jpj,jpk), &
874         & zrotn_adout(jpi,jpj,jpk),  &
875         & znu(jpi,jpj,jpk),          &
876         & znv(jpi,jpj,jpk)           &
877         & )
878
879
880      !==================================================================
881      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
882      !    dy = ( hdivb_tl, hdivn_tl )
883      !==================================================================
884
885      !--------------------------------------------------------------------
886      ! Reset the tangent and adjoint variables
887      !--------------------------------------------------------------------
888
889      zun_tlin    (:,:,:) = 0.0_wp
890      zvn_tlin    (:,:,:) = 0.0_wp
891      zhdivn_tlin (:,:,:) = 0.0_wp
892      zrotn_tlin  (:,:,:) = 0.0_wp
893      zhdivb_tlout(:,:,:) = 0.0_wp
894      zhdivn_tlout(:,:,:) = 0.0_wp
895      zrotb_tlout (:,:,:) = 0.0_wp
896      zrotn_tlout (:,:,:) = 0.0_wp
897      zhdivb_adin (:,:,:) = 0.0_wp
898      zhdivn_adin (:,:,:) = 0.0_wp
899      zrotb_adin  (:,:,:) = 0.0_wp
900      zrotn_adin  (:,:,:) = 0.0_wp
901      zrotn_adout (:,:,:) = 0.0_wp
902      zhdivn_adout(:,:,:) = 0.0_wp
903      zun_adout   (:,:,:) = 0.0_wp
904      zvn_adout   (:,:,:) = 0.0_wp
905
906      un_tl   (:,:,:) = 0.0_wp
907      vn_tl   (:,:,:) = 0.0_wp
908      hdivb_tl(:,:,:) = 0.0_wp
909      hdivn_tl(:,:,:) = 0.0_wp
910      rotb_tl (:,:,:) = 0.0_wp
911      rotn_tl (:,:,:) = 0.0_wp
912      hdivb_ad(:,:,:) = 0.0_wp
913      hdivn_ad(:,:,:) = 0.0_wp
914      rotb_ad (:,:,:) = 0.0_wp
915      rotn_ad (:,:,:) = 0.0_wp
916      un_ad   (:,:,:) = 0.0_wp
917      vn_ad   (:,:,:) = 0.0_wp
918
919      !--------------------------------------------------------------------
920      ! Initialize the tangent input with random noise: dx
921      !--------------------------------------------------------------------
922
923      CALL grid_random( znu, 'U', 0.0_wp, stdu )
924
925      CALL grid_random( znv, 'V', 0.0_wp, stdv )
926
927      DO jk = 1, jpk
928         DO jj = nldj, nlej
929            DO ji = nldi, nlei
930               zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
931               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
932            END DO
933         END DO
934      END DO
935
936      un_tl(:,:,:) = zun_tlin(:,:,:)
937      vn_tl(:,:,:) = zvn_tlin(:,:,:)
938
939      CALL div_cur_tan( nit000 )  ! Generate noise for before hdiv/rot fields
940
941      DO jk = 1, jpk
942         DO jj = nldj, nlej
943            DO ji = nldi, nlei
944               zhdivn_tlin(ji,jj,jk) = 0.5_wp * hdivn_tl(ji,jj,jk)
945               zrotn_tlin (ji,jj,jk) = 0.5_wp * rotn_tl (ji,jj,jk)
946            END DO
947         END DO
948      END DO
949
950      un_tl   (:,:,:) = 0.0_wp
951      vn_tl   (:,:,:) = 0.0_wp
952      hdivb_tl(:,:,:) = 0.0_wp
953      hdivn_tl(:,:,:) = 0.0_wp
954      rotb_tl (:,:,:) = 0.0_wp
955      rotn_tl (:,:,:) = 0.0_wp
956
957      !--------------------------------------------------------------------
958      ! Call the tangent routine: dy = L dx
959      !--------------------------------------------------------------------
960
961      un_tl   (:,:,:) = zun_tlin   (:,:,:)
962      vn_tl   (:,:,:) = zvn_tlin   (:,:,:)
963      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
964
965      CALL div_cur_tan( nit000 )
966
967      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
968      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
969
970      !--------------------------------------------------------------------
971      ! Initialize the adjoint variables: dy^* = W dy
972      !--------------------------------------------------------------------
973      DO jk = 1, jpk
974        DO jj = nldj, nlej
975           DO ji = nldi, nlei
976              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
977                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
978                 &               * tmask(ji,jj,jk)
979              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
980                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
981                 &               * tmask(ji,jj,jk)
982            END DO
983         END DO
984      END DO
985
986      !--------------------------------------------------------------------
987      ! Compute the scalar product: ( L dx )^T W dy
988      !--------------------------------------------------------------------
989
990      zsp1_1 = DOT_PRODUCT( zhdivb_tlout, zhdivb_adin )
991      zsp1_2 = DOT_PRODUCT( zhdivn_tlout, zhdivn_adin )
992      zsp1   = zsp1_1 + zsp1_2
993
994      !--------------------------------------------------------------------
995      ! Call the adjoint routine: dx^* = L^T dy^*
996      !--------------------------------------------------------------------
997
998      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
999      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
1000      rotb_ad (:,:,:) = 0.0_wp
1001      rotn_ad (:,:,:) = 0.0_wp
1002
1003      CALL div_cur_adj( nit000 )
1004
1005      zun_adout   (:,:,:) = un_ad   (:,:,:)
1006      zvn_adout   (:,:,:) = vn_ad   (:,:,:)
1007      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
1008
1009      !--------------------------------------------------------------------
1010      ! Compute the scalar product: dx^T L^T W dy
1011      !--------------------------------------------------------------------
1012
1013      zsp2_1 = DOT_PRODUCT( zun_tlin,    zun_adout    )
1014      zsp2_2 = DOT_PRODUCT( zvn_tlin,    zvn_adout    )
1015      zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
1016      zsp2   = zsp2_1 + zsp2_2 + zsp2_3
1017
1018      cl_name = 'div_cur_adj T1'
1019      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1020
1021      !=============================================================
1022      ! 2) dx = ( un_tl, vn_tl, rotn_tl ) and
1023      !    dy = ( rotb_tl, rotn_tl )
1024      !=============================================================
1025
1026      !--------------------------------------------------------------------
1027      ! Reset the tangent and adjoint variables
1028      !--------------------------------------------------------------------
1029
1030      un_tl   (:,:,:) = 0.0_wp
1031      vn_tl   (:,:,:) = 0.0_wp
1032      hdivb_tl(:,:,:) = 0.0_wp
1033      hdivn_tl(:,:,:) = 0.0_wp
1034      rotb_tl (:,:,:) = 0.0_wp
1035      rotn_tl (:,:,:) = 0.0_wp
1036      hdivb_ad(:,:,:) = 0.0_wp
1037      hdivn_ad(:,:,:) = 0.0_wp
1038      rotb_ad (:,:,:) = 0.0_wp
1039      rotn_ad (:,:,:) = 0.0_wp
1040      un_ad   (:,:,:) = 0.0_wp
1041      vn_ad   (:,:,:) = 0.0_wp
1042
1043      !--------------------------------------------------------------------
1044      ! Call the tangent routine: dy = L dx
1045      !--------------------------------------------------------------------
1046
1047      un_tl  (:,:,:) = zun_tlin  (:,:,:)
1048      vn_tl  (:,:,:) = zvn_tlin  (:,:,:)
1049      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1050
1051      CALL div_cur_tan( nit000 )
1052
1053      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
1054      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
1055
1056      !--------------------------------------------------------------------
1057      ! Initialize the adjoint variables: dy^* = W dy
1058      !--------------------------------------------------------------------
1059
1060      DO jk = 1, jpk
1061        DO jj = nldj, nlej
1062           DO ji = nldi, nlei
1063              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
1064                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1065              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
1066                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1067            END DO
1068         END DO
1069      END DO
1070
1071      !--------------------------------------------------------------------
1072      ! Compute the scalar product: ( L dx )^T W dy
1073      !--------------------------------------------------------------------
1074
1075      zsp1_1 = DOT_PRODUCT( zrotb_tlout, zrotb_adin )
1076      zsp1_2 = DOT_PRODUCT( zrotn_tlout, zrotn_adin )
1077      zsp1   = zsp1_1 + zsp1_2
1078
1079      !--------------------------------------------------------------------
1080      ! Call the adjoint routine: dx^* = L^T dy^*
1081      !--------------------------------------------------------------------
1082
1083      rotb_ad (:,:,:) = zrotb_adin(:,:,:)
1084      rotn_ad (:,:,:) = zrotn_adin(:,:,:)
1085      hdivb_ad(:,:,:) = 0.0_wp
1086      hdivn_ad(:,:,:) = 0.0_wp
1087
1088      CALL div_cur_adj( nit000 )
1089
1090      zun_adout  (:,:,:) = un_ad  (:,:,:)
1091      zvn_adout  (:,:,:) = vn_ad  (:,:,:)
1092      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1093
1094      !--------------------------------------------------------------------
1095      ! Compute the scalar product: dx^T L^T W dy
1096      !--------------------------------------------------------------------
1097
1098      zsp2_1 = DOT_PRODUCT( zun_tlin,   zun_adout   )
1099      zsp2_2 = DOT_PRODUCT( zvn_tlin,   zvn_adout   )
1100      zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1101      zsp2   = zsp2_1 + zsp2_2 + zsp2_3
1102
1103      cl_name = 'div_cur_adj T2'
1104      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1105
1106      !=============================================================
1107      ! 3) dx = ( un_tl, vn_tl, rotn_tl, hdin_tl ) and
1108      !    dy = ( rotb_tl, rotn_tl, hdivn_tl, hdivb_tl )
1109      !=============================================================
1110
1111      !--------------------------------------------------------------------
1112      ! Reset the tangent and adjoint variables
1113      !--------------------------------------------------------------------
1114
1115      un_tl   (:,:,:) = 0.0_wp
1116      vn_tl   (:,:,:) = 0.0_wp
1117      hdivb_tl(:,:,:) = 0.0_wp
1118      hdivn_tl(:,:,:) = 0.0_wp
1119      rotb_tl (:,:,:) = 0.0_wp
1120      rotn_tl (:,:,:) = 0.0_wp
1121      hdivb_ad(:,:,:) = 0.0_wp
1122      hdivn_ad(:,:,:) = 0.0_wp
1123      rotb_ad (:,:,:) = 0.0_wp
1124      rotn_ad (:,:,:) = 0.0_wp
1125      un_ad   (:,:,:) = 0.0_wp
1126      vn_ad   (:,:,:) = 0.0_wp
1127
1128      !--------------------------------------------------------------------
1129      ! Call the tangent routine: dy = L dx
1130      !--------------------------------------------------------------------
1131
1132      un_tl  (:,:,:) = zun_tlin  (:,:,:)
1133      vn_tl  (:,:,:) = zvn_tlin  (:,:,:)
1134      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1135      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
1136
1137      CALL div_cur_tan( nit000 )
1138
1139      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
1140      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
1141      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
1142      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
1143
1144      !--------------------------------------------------------------------
1145      ! Initialize the adjoint variables: dy^* = W dy
1146      !--------------------------------------------------------------------
1147      DO jk = 1, jpk
1148        DO jj = nldj, nlej
1149           DO ji = nldi, nlei
1150              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
1151                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1152                 &               * tmask(ji,jj,jk)
1153              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
1154                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1155                 &               * tmask(ji,jj,jk)
1156            END DO
1157         END DO
1158      END DO
1159      DO jk = 1, jpk
1160        DO jj = nldj, nlej
1161           DO ji = nldi, nlei
1162              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
1163                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1164              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
1165                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1166            END DO
1167         END DO
1168      END DO
1169
1170      !--------------------------------------------------------------------
1171      ! Compute the scalar product: ( L dx )^T W dy
1172      !--------------------------------------------------------------------
1173
1174      zsp1_1 = DOT_PRODUCT( zhdivb_tlout, zhdivb_adin )
1175      zsp1_2 = DOT_PRODUCT( zhdivn_tlout, zhdivn_adin )
1176      zsp1_3 = DOT_PRODUCT( zrotb_tlout, zrotb_adin )
1177      zsp1_4 = DOT_PRODUCT( zrotn_tlout, zrotn_adin )
1178      zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4
1179
1180      !--------------------------------------------------------------------
1181      ! Call the adjoint routine: dx^* = L^T dy^*
1182      !--------------------------------------------------------------------
1183
1184      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
1185      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
1186      rotb_ad (:,:,:) = zrotb_adin(:,:,:)
1187      rotn_ad (:,:,:) = zrotn_adin(:,:,:)
1188
1189      CALL div_cur_adj( nit000 )
1190
1191      zun_adout  (:,:,:) = un_ad  (:,:,:)
1192      zvn_adout  (:,:,:) = vn_ad  (:,:,:)
1193      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1194      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
1195
1196      !--------------------------------------------------------------------
1197      ! Compute the scalar product: dx^T L^T W dy
1198      !--------------------------------------------------------------------
1199
1200      zsp2_1 = DOT_PRODUCT( zun_tlin,   zun_adout   )
1201      zsp2_2 = DOT_PRODUCT( zvn_tlin,   zvn_adout   )
1202      zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1203      zsp2_4 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
1204      zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
1205
1206      cl_name = 'div_cur_adj T3'
1207      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1208
1209
1210      DEALLOCATE( &
1211         & zun_tlin,     &
1212         & zvn_tlin,     &
1213         & zhdivn_tlin,  &
1214         & zrotn_tlin,   &
1215         & zhdivb_tlout, &
1216         & zhdivn_tlout, &
1217         & zrotb_tlout,  &
1218         & zrotn_tlout,  &
1219         & zhdivb_adin,  &
1220         & zhdivn_adin,  &
1221         & zrotb_adin,   &
1222         & zrotn_adin,   &
1223         & zun_adout,    &
1224         & zvn_adout,    &
1225         & zhdivn_adout, &
1226         & zrotn_adout,  &
1227         & znu,          &
1228         & znv           &
1229         & )
1230
1231   END SUBROUTINE div_cur_adj_tst
1232#endif
1233
1234   !!======================================================================
1235
1236END MODULE divcur_tam
Note: See TracBrowser for help on using the repository browser.