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

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/divcur_tam.F90 @ 3627

Last change on this file since 3627 was 3627, checked in by rblod, 11 years ago

Correct long lines in TAM 4.3 see ticket #1010
  • Property svn:executable set to *
File size: 53.1 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      !!----------------------------------------------------------------------
358      IF( nn_timing == 1 )  CALL timing_start('div_cur_adj')
359      !
360      CALL wrk_alloc( jpi  , jpj+2, zwu               )
361      CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 )
362      !
363      IF( kt == nitend ) THEN
364         IF(lwp) WRITE(numout,*)
365         IF(lwp) WRITE(numout,*) 'div_cur_adj : horizontal velocity', &
366            &                    ' divergence and relative vorticity'
367         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   NOT optimal for auto-tasking case'
368      ENDIF
369
370      ! 4. Lateral boundary conditions on hdivn and rotn
371      ! ---------------------------------=======---======
372      CALL lbc_lnk_adj( rotn_ad , 'F', 1.0_wp )   ! F-point, no sign change
373      CALL lbc_lnk_adj( hdivn_ad, 'T', 1.0_wp )   ! T-point, no sign change
374
375      IF( nn_cla == 1 )   CALL cla_div_adj    ( kt )             ! Cross Land Advection (Update Hor. divergence)
376      IF( ln_rnf      )   CALL sbc_rnf_div_adj( hdivn_ad )       ! runoffs (update hdivn field)
377      !                                                ! ===============
378      DO jk = jpkm1, 1, -1                             ! Horizontal slab
379         !                                             ! ===============
380         ! local adjoint workspace initialization
381         zwu(:,:) = 0.0_wp
382         zwv(:,:) = 0.0_wp
383         !                                             ! --------
384         ! relative vorticity                          !   rot
385         !                                             ! --------
386         DO jl = npcoa(4,jk), 1, -1
387            ii = nicoa(jl,4,jk)
388            ij = njcoa(jl,4,jk)
389            rotn_ad(ii,ij,jk) = -1.0_wp * rotn_ad(ii,ij,jk) &
390               &                / ( e1f(ii,ij) * e2f(ii,ij) )
391            zwu(ii,ij  ) = zwu(ii,ij  ) - 4.0_wp * rotn_ad(ii,ij,jk)
392            zwu(ii,ij-1) = zwu(ii,ij-1) +          rotn_ad(ii,ij,jk)
393            zwu(ii,ij-2) = zwu(ii,ij-2) - 0.2_wp * rotn_ad(ii,ij,jk)
394            rotn_ad(ii,ij,jk) = 0.0_wp
395         END DO
396         DO jl = npcoa(3,jk), 1, -1
397            ii = nicoa(jl,3,jk)
398            ij = njcoa(jl,3,jk)
399            rotn_ad(ii,ij,jk) = -1.0_wp * rotn_ad(ii,ij,jk) &
400               &                / ( e1f(ii,ij) * e2f(ii,ij) )
401            zwu(ii,ij+1) = zwu(ii,ij+1) + 4.0_wp * rotn_ad(ii,ij,jk)
402            zwu(ii,ij+2) = zwu(ii,ij+2) -          rotn_ad(ii,ij,jk)
403            zwu(ii,ij+3) = zwu(ii,ij+3) + 0.2_wp * rotn_ad(ii,ij,jk)
404            rotn_ad(ii,ij,jk) = 0.0_wp
405         END DO
406         DO jl = npcoa(2,jk), 1, -1
407            ii = nicoa(jl,2,jk)
408            ij = njcoa(jl,2,jk)
409            rotn_ad(ii,ij,jk) =  1.0_wp * rotn_ad(ii,ij,jk) &
410               &                / ( e1f(ii,ij) * e2f(ii,ij) )
411            zwv(ii  ,ij) = zwv(ii  ,ij) - 4.0_wp * rotn_ad(ii,ij,jk)
412            zwv(ii-1,ij) = zwv(ii-1,ij) +          rotn_ad(ii,ij,jk)
413            zwv(ii-2,ij) = zwv(ii-2,ij) - 0.2_wp * rotn_ad(ii,ij,jk)
414            rotn_ad(ii,ij,jk) = 0.0_wp
415         END DO
416         ! second order accurate scheme along straight coast
417         DO jl = npcoa(1,jk), 1, -1
418            ii = nicoa(jl,1,jk)
419            ij = njcoa(jl,1,jk)
420            rotn_ad(ii,ij,jk) =  1.0_wp * rotn_ad(ii,ij,jk) &
421               &                / ( e1f(ii,ij) * e2f(ii,ij) )
422            zwv(ii+1,ij) = zwv(ii+1,ij) + 4.0_wp * rotn_ad(ii,ij,jk)
423            zwv(ii+2,ij) = zwv(ii+2,ij) -          rotn_ad(ii,ij,jk)
424            zwv(ii+3,ij) = zwv(ii+3,ij) + 0.2_wp * rotn_ad(ii,ij,jk)
425            rotn_ad(ii,ij,jk) = 0.0_wp
426         END DO
427         ! relative vorticity (vertical component of the velocity curl)
428         DO jj = jpjm1, 1, -1
429            DO ji = fs_jpim1, 1, -1   ! vector opt.
430               rotn_ad(ji,jj,jk) = rotn_ad(ji,jj,jk) * fmask(ji,jj,jk) &
431                  &                / ( e1f(ji,jj) * e2f(ji,jj) )
432               zwv(ji  ,jj  ) = zwv(ji  ,jj  ) - rotn_ad(ji,jj,jk)
433               zwu(ji  ,jj  ) = zwu(ji  ,jj  ) + rotn_ad(ji,jj,jk)
434               zwu(ji  ,jj+1) = zwu(ji  ,jj+1) - rotn_ad(ji,jj,jk)
435               zwv(ji+1,jj  ) = zwv(ji+1,jj  ) + rotn_ad(ji,jj,jk)
436               rotn_ad(ji,jj,jk) = 0.0
437            END DO
438         END DO
439         ! North-South boundary conditions
440         IF( nperio == 3 .OR. nperio == 4 ) THEN
441            ! north fold ( Grid defined with a T-point pivot) ORCA 2 degree
442            DO ji = jpi-1, 1, -1
443               iju = jpi - ji + 1
444               zwu(iju,jpj-4) = zwu(iju,jpj-4) - zwu(ji,jpj+2)
445               zwu(ji ,jpj+2) = 0.0_wp
446               zwu(iju,jpj-3) = zwu(iju,jpj-3) - zwu(ji,jpj+1)
447               zwu(ji ,jpj+1) = 0.0_wp
448            END DO
449            zwu(jpi,jpj+2) = 0.0_wp
450            zwu(jpi,jpj+1) = 0.0_wp
451         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
452            ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degree
453            DO ji = jpi+2, jpi/2+1, -1
454               ijt = jpi - ji + 1
455               zwv(ijt,jpjm1) = zwv(ijt,jpjm1) - zwv(ji,jpjm1)
456               zwv(ji ,jpjm1) = 0.0_wp
457            END DO
458            DO ji = jpi+2, -1, -1
459               ijt = jpi - ji + 1
460               zwv(ijt,jpj-2) = zwv(ijt,jpj-2) - zwv(ji,jpj  )
461               zwv(ji ,jpj  ) = 0.0_wp
462            END DO
463            DO ji = jpi-1, 1, -1
464               iju = jpi - ji
465               zwu(iju,jpj-3) = zwu(iju,jpj-3) - zwu(ji,jpj+2)
466               zwu(ji ,jpj+2) = 0.0_wp
467               zwu(iju,jpj-2) = zwu(iju,jpj-2) - zwu(ji,jpj+1)
468               zwu(ji ,jpj+1) = 0.0_wp
469               zwu(iju,jpj-1) = zwu(iju,jpj-1) - zwu(ji,jpj  )
470               zwu(ji ,jpj  ) = 0.0_wp
471            END DO
472            zwu(jpi,jpj+2) = 0.0_wp
473            zwu(jpi,jpj+1) = 0.0_wp
474         ELSE
475            ! closed
476            zwu(:,jpj+2) = 0.0_wp
477            zwu(:,jpj+1) = 0.0_wp
478         ENDIF
479         ! East-West boundary conditions
480         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
481            zwv(  4  ,:) = zwv(  4  ,:) + zwv(jpi+2,:)
482            zwv(jpi+2,:) = 0.0_wp
483            zwv(  3  ,:) = zwv(  3  ,:) + zwv(jpi+1,:)
484            zwv(jpi+1,:) = 0.0_wp
485            zwv(jpi-3,:) = zwv(jpi-3,:) + zwv( -1  ,:)
486            zwv( -1  ,:) = 0.0_wp
487            zwv(jpi-2,:) = zwv(jpi-2,:) + zwv(  0  ,:)
488            zwv(  0  ,:) = 0.0_wp
489         ELSE
490            zwv(jpi+2,:) = 0.0_wp
491            zwv(jpi+1,:) = 0.0_wp
492            zwv( -1  ,:) = 0.0_wp
493            zwv(  0  ,:) = 0.0_wp
494         ENDIF
495         ! contravariant velocity (extended for lateral b.c.)
496         ! inside the model domain
497         DO jj = jpj, 1, -1
498            DO ji = jpi, 1, -1
499               vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e2v(ji,jj) * zwv(ji,jj)
500               un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e1u(ji,jj) * zwu(ji,jj)
501            END DO
502         END DO
503         !
504         IF( .NOT. AGRIF_Root() ) THEN
505            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn_ad(nlci-1 , :     ,jk) = 0.0_wp      ! east
506            IF ((nbondi == -1).OR.(nbondi == 2)) hdivn_ad(2      , :     ,jk) = 0.0_wp      ! west
507            IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn_ad(:      ,nlcj-1 ,jk) = 0.0_wp      ! north
508            IF ((nbondj == -1).OR.(nbondj == 2)) hdivn_ad(:      ,2      ,jk) = 0.0_wp      ! south
509         ENDIF
510         !                                             ! --------
511         ! Horizontal divergence                       !   div
512         !                                             ! --------
513         DO jj = jpjm1, 2, -1
514            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
515               hdivn_ad(ji,jj,jk) =  hdivn_ad(ji,jj,jk)  &
516                  &                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
517               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
518                  &                  + e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) &
519                  &                                   * hdivn_ad(ji,jj,jk)
520               un_ad(ji-1,jj  ,jk) = un_ad(ji-1,jj  ,jk) &
521                  &                  - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) &
522                  &                                   * hdivn_ad(ji,jj,jk)
523               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
524                  &                  + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) &
525                  &                                   * hdivn_ad(ji,jj,jk)
526               vn_ad(ji  ,jj-1,jk) = vn_ad(ji  ,jj-1,jk) &
527                  &                  - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) &
528                  &                                   * hdivn_ad(ji,jj,jk)
529               hdivn_ad(ji,jj,jk) = 0.0_wp
530            END DO
531         END DO
532         rotn_ad (:,:,jk) = rotn_ad (:,:,jk) + rotb_ad (:,:,jk)  ! time swap
533         rotb_ad (:,:,jk) = 0.0_wp
534         hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + hdivb_ad(:,:,jk)  ! time swap
535         hdivb_ad(:,:,jk) = 0.0_wp
536         !                                             ! ===============
537      END DO                                           !   End of slab
538      !                                                ! ===============
539      CALL wrk_dealloc( jpi  , jpj+2, zwu               )
540      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 )
541      !
542      IF( nn_timing == 1 )  CALL timing_stop('div_cur_adj')
543
544   END SUBROUTINE div_cur_adj
545
546#else
547   !!----------------------------------------------------------------------
548   !!   Default option                           2nd order centered schemes
549   !!----------------------------------------------------------------------
550   SUBROUTINE div_cur_tan( kt )
551      !!----------------------------------------------------------------------
552      !!                  ***  ROUTINE div_cur_tan  ***
553      !!
554      !! ** Purpose of direct routine :
555      !!      compute the horizontal divergence and the relative
556      !!      vorticity at before and now time-step
557      !!
558      !! ** Method of direct routine :
559      !!              - Divergence:
560      !!      - save the divergence computed at the previous time-step
561      !!      (note that the Asselin filter has not been applied on hdivb)
562      !!      - compute the now divergence given by :
563      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
564      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
565      !!      above expression
566      !!      - apply lateral boundary conditions on hdivn
567      !!              - Relavtive Vorticity :
568      !!      - save the curl computed at the previous time-step (rotb = rotn)
569      !!      (note that the Asselin time filter has not been applied to rotb)
570      !!      - compute the now curl in tensorial formalism:
571      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
572      !!      - apply lateral boundary conditions on rotn through a call to
573      !!      routine lbc_lnk routine.
574      !!      Note: Coastal boundary condition: lateral friction set through
575      !!      the value of fmask along the coast (see dommsk.F90) and shlat
576      !!      (namelist parameter)
577      !!
578      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
579      !!              - update rotb , rotn , the before & now rel. vorticity
580      !!
581      !! History of the direct routine:
582      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
583      !!   4.0  !  91-11  (G. Madec)
584      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
585      !!   7.0  !  96-01  (G. Madec)  s-coordinates
586      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
587      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
588      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
589      !!        !  05-01  (J. Chanut) Unstructured open boundaries
590      !! History of the TAM routine:
591      !!   9.0  !  08-06  (A. Vidard) Skeleton
592      !!        !  08-07  (A. Weaver) tangent of the 02-09 version
593      !!        !  08-11  (A. Vidard) tangent of the 05-01 version
594      !!----------------------------------------------------------------------
595      !! * Arguments
596      INTEGER, INTENT( in ) :: &
597         & kt     ! ocean time-step index
598
599      !! * Local declarations
600      INTEGER :: &
601         & ji,  & ! dummy loop indices
602         & jj,  &
603         & jk
604      !!----------------------------------------------------------------------
605      IF( nn_timing == 1 )  CALL timing_start('div_cur_tan')
606      !
607      IF( kt == nit000 ) THEN
608         IF(lwp) WRITE(numout,*)
609         IF(lwp) WRITE(numout,*) 'div_cur_tan : horizontal velocity divergence and'
610         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   relative vorticity'
611      ENDIF
612      !                                                ! ===============
613      DO jk = 1, jpkm1                                 ! Horizontal slab
614         !                                             ! ===============
615         hdivb_tl(:,:,jk) = hdivn_tl(:,:,jk)    ! time swap of div arrays
616         rotb_tl (:,:,jk) = rotn_tl (:,:,jk)    ! time swap of rot arrays
617         !                                             ! --------
618         ! Horizontal divergence                       !   div
619         !                                             ! --------
620         DO jj = 2, jpjm1
621            DO ji = fs_2, fs_jpim1   ! vector opt.
622               hdivn_tl(ji,jj,jk) =   &
623                  & (   e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un_tl(ji  ,jj  ,jk) &
624                  &   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un_tl(ji-1,jj  ,jk) &
625                  &   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn_tl(ji  ,jj  ,jk) &
626                  &   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn_tl(ji  ,jj-1,jk) &
627                  & ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
628            END DO
629         END DO
630         !                                             ! --------
631         ! relative vorticity                          !   rot
632         !                                             ! --------
633         DO jj = 1, jpjm1
634            DO ji = 1, fs_jpim1   ! vector opt.
635               rotn_tl(ji,jj,jk) = (   e2v(ji+1,jj  ) * vn_tl(ji+1,jj  ,jk) &
636                  &                  - e2v(ji  ,jj  ) * vn_tl(ji  ,jj  ,jk) &
637                  &                  - e1u(ji  ,jj+1) * un_tl(ji  ,jj+1,jk) &
638                  &                  + e1u(ji  ,jj  ) * un_tl(ji  ,jj  ,jk) &
639                  &                ) * fmask(ji,jj,jk) / ( e1f(ji,jj) * e2f(ji,jj) )
640            END DO
641         END DO
642         !                                             ! ===============
643      END DO                                           !   End of slab
644      !                                                ! ===============
645      IF( ln_rnf      )   CALL sbc_rnf_div_tan( hdivn_tl )       ! runoffs (update hdivn field)
646      IF( nn_cla == 1 )   CALL cla_div_tan    ( kt )             ! Cross Land Advection (update hdivn field)
647      !!
648      CALL lbc_lnk( hdivn_tl, 'T', 1. )
649      CALL lbc_lnk( rotn_tl , 'F', 1. )     ! lateral boundary cond. (no sign change)
650      !
651      IF( nn_timing == 1 )  CALL timing_stop('div_cur_tan')
652   END SUBROUTINE div_cur_tan
653
654   SUBROUTINE div_cur_adj( kt )
655      !!----------------------------------------------------------------------
656      !!                  ***  ROUTINE div_cur_adj  ***
657      !!
658      !! ** Purpose of direct routine :
659      !!      compute the horizontal divergence and the relative
660      !!      vorticity at before and now time-step
661      !!
662      !! ** Method of direct routine :
663      !!              - Divergence:
664      !!      - save the divergence computed at the previous time-step
665      !!      (note that the Asselin filter has not been applied on hdivb)
666      !!      - compute the now divergence given by :
667      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
668      !!      Note: if lk_zco=T, e3u=e3v=e3t, they are simplified in the
669      !!      above expression
670      !!      - apply lateral boundary conditions on hdivn
671      !!              - Relavtive Vorticity :
672      !!      - save the curl computed at the previous time-step (rotb = rotn)
673      !!      (note that the Asselin time filter has not been applied to rotb)
674      !!      - compute the now curl in tensorial formalism:
675      !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
676      !!      - apply lateral boundary conditions on rotn through a call to
677      !!      routine lbc_lnk routine.
678      !!      Note: Coastal boundary condition: lateral friction set through
679      !!      the value of fmask along the coast (see dommsk.F90) and shlat
680      !!      (namelist parameter)
681      !!
682      !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
683      !!              - update rotb , rotn , the before & now rel. vorticity
684      !!
685      !! History of the direct routine:
686      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
687      !!   4.0  !  91-11  (G. Madec)
688      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
689      !!   7.0  !  96-01  (G. Madec)  s-coordinates
690      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
691      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
692      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
693      !! History of the TAM routine:
694      !!   9.0  !  08-06  (A. Vidard) Skeleton
695      !!   9.0  !  08-07  (A. Weaver)
696      !!----------------------------------------------------------------------
697      !! * Arguments
698      INTEGER, INTENT( in ) :: &
699         & kt     ! ocean time-step index
700
701      !! * Local declarations
702      INTEGER :: &
703         & ji,  & ! dummy loop indices
704         & jj,  &
705         & jk
706      !!----------------------------------------------------------------------
707      !
708      if( nn_timing == 1 )  call timing_start('div_cur_adj')
709      !
710      IF( kt == nitend ) THEN
711         IF(lwp) WRITE(numout,*)
712         IF(lwp) WRITE(numout,*) 'div_cur_adj : horizontal velocity divergence and'
713         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   relative vorticity'
714      ENDIF
715      ! 4. Lateral boundary conditions on hdivn and rotn
716      ! ---------------------------------=======---======
717      CALL lbc_lnk_adj( rotn_ad , 'F', 1.0_wp )     ! F-point, no sign change
718      CALL lbc_lnk_adj( hdivn_ad, 'T', 1.0_wp )     ! T-point, no sign change
719      !!
720      IF( nn_cla == 1 )   CALL cla_div_adj    ( kt )             ! Cross Land Advection (update hdivn field)
721      IF( ln_rnf      )   CALL sbc_rnf_div_adj( hdivn_ad )       ! runoffs (update hdivn field)
722      !                                                ! ===============
723      DO jk = jpkm1, 1, -1                             ! Horizontal slab
724         !                                             ! ===============
725         !                                             ! --------
726         ! relative vorticity                          !   rot
727         !                                             ! --------
728         DO jj = jpjm1, 1, -1
729            DO ji = fs_jpim1, 1, -1   ! vector opt.
730               rotn_ad(ji,jj,jk) =  rotn_ad(ji,jj,jk) * fmask(ji,jj,jk) &
731                  &                       / ( e1f(ji,jj) * e2f(ji,jj) )
732               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
733                  &                  + e1u(ji  ,jj  ) * rotn_ad(ji,jj,jk)
734               un_ad(ji  ,jj+1,jk) = un_ad(ji  ,jj+1,jk) &
735                  &                  - e1u(ji  ,jj+1) * rotn_ad(ji,jj,jk)
736               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
737                  &                  - e2v(ji  ,jj  ) * rotn_ad(ji,jj,jk)
738               vn_ad(ji+1,jj  ,jk) = vn_ad(ji+1,jj  ,jk) &
739                  &                  + e2v(ji+1,jj  ) * rotn_ad(ji,jj,jk)
740               rotn_ad(ji,jj,jk) = 0.0_wp
741            END DO
742         END DO
743         !                                             ! --------
744         ! Horizontal divergence                       !   div
745         !                                             ! --------
746         DO jj = jpjm1, 2, -1
747            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
748               hdivn_ad(ji,jj,jk) =  hdivn_ad(ji,jj,jk)  &
749                  &                  / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
750               un_ad(ji  ,jj  ,jk) = un_ad(ji  ,jj  ,jk) &
751                  &                  + e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) &
752                  &                                   * hdivn_ad(ji,jj,jk)
753               un_ad(ji-1,jj  ,jk) = un_ad(ji-1,jj  ,jk) &
754                  &                  - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) &
755                  &                                   * hdivn_ad(ji,jj,jk)
756               vn_ad(ji  ,jj  ,jk) = vn_ad(ji  ,jj  ,jk) &
757                  &                  + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) &
758                  &                                   * hdivn_ad(ji,jj,jk)
759               vn_ad(ji  ,jj-1,jk) = vn_ad(ji  ,jj-1,jk) &
760                  &                  - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) &
761                  &                                   * hdivn_ad(ji,jj,jk)
762               hdivn_ad(ji,jj,jk) = 0.0_wp
763            END DO
764         END DO
765         !
766         rotn_ad (:,:,jk) = rotn_ad (:,:,jk) + rotb_ad (:,:,jk)  ! time swap
767         rotb_ad (:,:,jk) = 0.0_wp
768         hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + hdivb_ad(:,:,jk)  ! time swap
769         hdivb_ad(:,:,jk) = 0.0_wp
770         !                                             ! ===============
771      END DO                                           !   End of slab
772      !                                                ! ===============
773      if( nn_timing == 1 )  call timing_stop('div_cur_adj')
774      !
775   END SUBROUTINE div_cur_adj
776
777#endif
778
779   SUBROUTINE div_cur_adj_tst( kumadt )
780      !!-----------------------------------------------------------------------
781      !!
782      !!          ***  ROUTINE div_cur_adj_tst : TEST OF div_cur_adj  ***
783      !!
784      !! ** Purpose : Test the adjoint routine.
785      !!
786      !! ** Method  : Verify the scalar product
787      !!
788      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
789      !!
790      !!              where  L   = tangent routine
791      !!                     L^T = adjoint routine
792      !!                     W   = diagonal matrix of scale factors
793      !!                     dx  = input perturbation (random field)
794      !!                     dy  = L dx
795      !!
796      !! ** Action  : Separate tests are applied for the following dx and dy:
797      !!
798      !!      1) dx = ( un_tl, vn_tl ) and
799      !!         dy = ( hdivn_tl )
800      !!      2) dx = ( un_tl, vn_tl ) and
801      !!         dy = ( rotntl )
802      !!
803      !! History :
804      !!        ! 08-07 (A. Weaver)
805      !!-----------------------------------------------------------------------
806
807      !! * Modules used
808      !! * Arguments
809      INTEGER, INTENT(IN) :: &
810         & kumadt             ! Output unit
811
812      INTEGER :: &
813         & ji,    &        ! dummy loop indices
814         & jj,    &
815         & jk
816
817      !! * Local declarations
818      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
819         & zun_tlin,     & ! Tangent input: now u-velocity
820         & zvn_tlin,     & ! Tangent input: now v-velocity
821         & zhdivn_tlin,  & ! Tangent input: now horizontal divergence
822         & zrotn_tlin,   & ! Tangent input: now relative vorticity
823         & zhdivb_tlout, & ! Tangent output: before horizontal divergence
824         & zhdivn_tlout, & ! Tangent output: now horizontal divergence
825         & zrotb_tlout,  & ! Tangent output: before relative vorticity
826         & zrotn_tlout,  & ! Tangent output: now relative vorticity
827         & zhdivb_adin,  & ! Adjoint input: before horizontal divergence
828         & zhdivn_adin,  & ! Adjoint input: now horizontal divergence
829         & zrotb_adin,   & ! Adjoint input: before relative vorticity
830         & zrotn_adin,   & ! Adjoint input: now relative vorticity
831         & zun_adout,    & ! Adjoint output: now u-velocity
832         & zvn_adout,    & ! Adjoint output: now v-velocity
833         & zhdivn_adout, & ! Adjoint output: now horizontal divergence
834         & zrotn_adout,  & ! Adjoint output: now relative vorticity
835         & znu,          & ! 3D random field for u
836         & znv             ! 3D random field for v
837
838      REAL(KIND=wp) :: &
839                           ! random field standard deviation for:
840         & zsp1,         & ! scalar product involving the tangent routine
841         & zsp1_1,       & !   scalar product components
842         & zsp1_2,       &
843         & zsp1_3,       & !
844         & zsp1_4,       &
845         & zsp2,         & ! scalar product involving the adjoint routine
846         & zsp2_1,       & !   scalar product components
847         & zsp2_2,       &
848         & zsp2_3,       &
849         & zsp2_4
850
851      CHARACTER(LEN=14) :: cl_name
852
853      ! Allocate memory
854
855      ALLOCATE( &
856         & zun_tlin(jpi,jpj,jpk),     &
857         & zvn_tlin(jpi,jpj,jpk),     &
858         & zhdivn_tlin(jpi,jpj,jpk),  &
859         & zrotn_tlin(jpi,jpj,jpk),   &
860         & zhdivb_tlout(jpi,jpj,jpk), &
861         & zhdivn_tlout(jpi,jpj,jpk), &
862         & zrotb_tlout(jpi,jpj,jpk),  &
863         & zrotn_tlout(jpi,jpj,jpk),  &
864         & zhdivb_adin(jpi,jpj,jpk),  &
865         & zhdivn_adin(jpi,jpj,jpk),  &
866         & zrotb_adin(jpi,jpj,jpk),   &
867         & zrotn_adin(jpi,jpj,jpk),   &
868         & zun_adout(jpi,jpj,jpk),    &
869         & zvn_adout(jpi,jpj,jpk),    &
870         & zhdivn_adout(jpi,jpj,jpk), &
871         & zrotn_adout(jpi,jpj,jpk),  &
872         & znu(jpi,jpj,jpk),          &
873         & znv(jpi,jpj,jpk)           &
874         & )
875
876
877      !==================================================================
878      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
879      !    dy = ( hdivb_tl, hdivn_tl )
880      !==================================================================
881
882      !--------------------------------------------------------------------
883      ! Reset the tangent and adjoint variables
884      !--------------------------------------------------------------------
885
886      zun_tlin    (:,:,:) = 0.0_wp
887      zvn_tlin    (:,:,:) = 0.0_wp
888      zhdivn_tlin (:,:,:) = 0.0_wp
889      zrotn_tlin  (:,:,:) = 0.0_wp
890      zhdivb_tlout(:,:,:) = 0.0_wp
891      zhdivn_tlout(:,:,:) = 0.0_wp
892      zrotb_tlout (:,:,:) = 0.0_wp
893      zrotn_tlout (:,:,:) = 0.0_wp
894      zhdivb_adin (:,:,:) = 0.0_wp
895      zhdivn_adin (:,:,:) = 0.0_wp
896      zrotb_adin  (:,:,:) = 0.0_wp
897      zrotn_adin  (:,:,:) = 0.0_wp
898      zrotn_adout (:,:,:) = 0.0_wp
899      zhdivn_adout(:,:,:) = 0.0_wp
900      zun_adout   (:,:,:) = 0.0_wp
901      zvn_adout   (:,:,:) = 0.0_wp
902
903      un_tl   (:,:,:) = 0.0_wp
904      vn_tl   (:,:,:) = 0.0_wp
905      hdivb_tl(:,:,:) = 0.0_wp
906      hdivn_tl(:,:,:) = 0.0_wp
907      rotb_tl (:,:,:) = 0.0_wp
908      rotn_tl (:,:,:) = 0.0_wp
909      hdivb_ad(:,:,:) = 0.0_wp
910      hdivn_ad(:,:,:) = 0.0_wp
911      rotb_ad (:,:,:) = 0.0_wp
912      rotn_ad (:,:,:) = 0.0_wp
913      un_ad   (:,:,:) = 0.0_wp
914      vn_ad   (:,:,:) = 0.0_wp
915
916      !--------------------------------------------------------------------
917      ! Initialize the tangent input with random noise: dx
918      !--------------------------------------------------------------------
919
920      CALL grid_random( znu, 'U', 0.0_wp, stdu )
921
922      CALL grid_random( znv, 'V', 0.0_wp, stdv )
923
924      DO jk = 1, jpk
925         DO jj = nldj, nlej
926            DO ji = nldi, nlei
927               zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
928               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
929            END DO
930         END DO
931      END DO
932
933      un_tl(:,:,:) = zun_tlin(:,:,:)
934      vn_tl(:,:,:) = zvn_tlin(:,:,:)
935
936      CALL div_cur_tan( nit000 )  ! Generate noise for before hdiv/rot fields
937
938      DO jk = 1, jpk
939         DO jj = nldj, nlej
940            DO ji = nldi, nlei
941               zhdivn_tlin(ji,jj,jk) = 0.5_wp * hdivn_tl(ji,jj,jk)
942               zrotn_tlin (ji,jj,jk) = 0.5_wp * rotn_tl (ji,jj,jk)
943            END DO
944         END DO
945      END DO
946
947      un_tl   (:,:,:) = 0.0_wp
948      vn_tl   (:,:,:) = 0.0_wp
949      hdivb_tl(:,:,:) = 0.0_wp
950      hdivn_tl(:,:,:) = 0.0_wp
951      rotb_tl (:,:,:) = 0.0_wp
952      rotn_tl (:,:,:) = 0.0_wp
953
954      !--------------------------------------------------------------------
955      ! Call the tangent routine: dy = L dx
956      !--------------------------------------------------------------------
957
958      un_tl   (:,:,:) = zun_tlin   (:,:,:)
959      vn_tl   (:,:,:) = zvn_tlin   (:,:,:)
960      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
961
962      CALL div_cur_tan( nit000 )
963
964      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
965      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
966
967      !--------------------------------------------------------------------
968      ! Initialize the adjoint variables: dy^* = W dy
969      !--------------------------------------------------------------------
970      DO jk = 1, jpk
971        DO jj = nldj, nlej
972           DO ji = nldi, nlei
973              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
974                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
975                 &               * tmask(ji,jj,jk)
976              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
977                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
978                 &               * tmask(ji,jj,jk)
979            END DO
980         END DO
981      END DO
982
983      !--------------------------------------------------------------------
984      ! Compute the scalar product: ( L dx )^T W dy
985      !--------------------------------------------------------------------
986
987      zsp1_1 = DOT_PRODUCT( zhdivb_tlout, zhdivb_adin )
988      zsp1_2 = DOT_PRODUCT( zhdivn_tlout, zhdivn_adin )
989      zsp1   = zsp1_1 + zsp1_2
990
991      !--------------------------------------------------------------------
992      ! Call the adjoint routine: dx^* = L^T dy^*
993      !--------------------------------------------------------------------
994
995      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
996      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
997      rotb_ad (:,:,:) = 0.0_wp
998      rotn_ad (:,:,:) = 0.0_wp
999
1000      CALL div_cur_adj( nit000 )
1001
1002      zun_adout   (:,:,:) = un_ad   (:,:,:)
1003      zvn_adout   (:,:,:) = vn_ad   (:,:,:)
1004      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
1005
1006      !--------------------------------------------------------------------
1007      ! Compute the scalar product: dx^T L^T W dy
1008      !--------------------------------------------------------------------
1009
1010      zsp2_1 = DOT_PRODUCT( zun_tlin,    zun_adout    )
1011      zsp2_2 = DOT_PRODUCT( zvn_tlin,    zvn_adout    )
1012      zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
1013      zsp2   = zsp2_1 + zsp2_2 + zsp2_3
1014
1015      cl_name = 'div_cur_adj T1'
1016      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1017
1018      !=============================================================
1019      ! 2) dx = ( un_tl, vn_tl, rotn_tl ) and
1020      !    dy = ( rotb_tl, rotn_tl )
1021      !=============================================================
1022
1023      !--------------------------------------------------------------------
1024      ! Reset the tangent and adjoint variables
1025      !--------------------------------------------------------------------
1026
1027      un_tl   (:,:,:) = 0.0_wp
1028      vn_tl   (:,:,:) = 0.0_wp
1029      hdivb_tl(:,:,:) = 0.0_wp
1030      hdivn_tl(:,:,:) = 0.0_wp
1031      rotb_tl (:,:,:) = 0.0_wp
1032      rotn_tl (:,:,:) = 0.0_wp
1033      hdivb_ad(:,:,:) = 0.0_wp
1034      hdivn_ad(:,:,:) = 0.0_wp
1035      rotb_ad (:,:,:) = 0.0_wp
1036      rotn_ad (:,:,:) = 0.0_wp
1037      un_ad   (:,:,:) = 0.0_wp
1038      vn_ad   (:,:,:) = 0.0_wp
1039
1040      !--------------------------------------------------------------------
1041      ! Call the tangent routine: dy = L dx
1042      !--------------------------------------------------------------------
1043
1044      un_tl  (:,:,:) = zun_tlin  (:,:,:)
1045      vn_tl  (:,:,:) = zvn_tlin  (:,:,:)
1046      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1047
1048      CALL div_cur_tan( nit000 )
1049
1050      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
1051      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
1052
1053      !--------------------------------------------------------------------
1054      ! Initialize the adjoint variables: dy^* = W dy
1055      !--------------------------------------------------------------------
1056
1057      DO jk = 1, jpk
1058        DO jj = nldj, nlej
1059           DO ji = nldi, nlei
1060              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
1061                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1062              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
1063                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1064            END DO
1065         END DO
1066      END DO
1067
1068      !--------------------------------------------------------------------
1069      ! Compute the scalar product: ( L dx )^T W dy
1070      !--------------------------------------------------------------------
1071
1072      zsp1_1 = DOT_PRODUCT( zrotb_tlout, zrotb_adin )
1073      zsp1_2 = DOT_PRODUCT( zrotn_tlout, zrotn_adin )
1074      zsp1   = zsp1_1 + zsp1_2
1075
1076      !--------------------------------------------------------------------
1077      ! Call the adjoint routine: dx^* = L^T dy^*
1078      !--------------------------------------------------------------------
1079
1080      rotb_ad (:,:,:) = zrotb_adin(:,:,:)
1081      rotn_ad (:,:,:) = zrotn_adin(:,:,:)
1082      hdivb_ad(:,:,:) = 0.0_wp
1083      hdivn_ad(:,:,:) = 0.0_wp
1084
1085      CALL div_cur_adj( nit000 )
1086
1087      zun_adout  (:,:,:) = un_ad  (:,:,:)
1088      zvn_adout  (:,:,:) = vn_ad  (:,:,:)
1089      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1090
1091      !--------------------------------------------------------------------
1092      ! Compute the scalar product: dx^T L^T W dy
1093      !--------------------------------------------------------------------
1094
1095      zsp2_1 = DOT_PRODUCT( zun_tlin,   zun_adout   )
1096      zsp2_2 = DOT_PRODUCT( zvn_tlin,   zvn_adout   )
1097      zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1098      zsp2   = zsp2_1 + zsp2_2 + zsp2_3
1099
1100      cl_name = 'div_cur_adj T2'
1101      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1102
1103      !=============================================================
1104      ! 3) dx = ( un_tl, vn_tl, rotn_tl, hdin_tl ) and
1105      !    dy = ( rotb_tl, rotn_tl, hdivn_tl, hdivb_tl )
1106      !=============================================================
1107
1108      !--------------------------------------------------------------------
1109      ! Reset the tangent and adjoint variables
1110      !--------------------------------------------------------------------
1111
1112      un_tl   (:,:,:) = 0.0_wp
1113      vn_tl   (:,:,:) = 0.0_wp
1114      hdivb_tl(:,:,:) = 0.0_wp
1115      hdivn_tl(:,:,:) = 0.0_wp
1116      rotb_tl (:,:,:) = 0.0_wp
1117      rotn_tl (:,:,:) = 0.0_wp
1118      hdivb_ad(:,:,:) = 0.0_wp
1119      hdivn_ad(:,:,:) = 0.0_wp
1120      rotb_ad (:,:,:) = 0.0_wp
1121      rotn_ad (:,:,:) = 0.0_wp
1122      un_ad   (:,:,:) = 0.0_wp
1123      vn_ad   (:,:,:) = 0.0_wp
1124
1125      !--------------------------------------------------------------------
1126      ! Call the tangent routine: dy = L dx
1127      !--------------------------------------------------------------------
1128
1129      un_tl  (:,:,:) = zun_tlin  (:,:,:)
1130      vn_tl  (:,:,:) = zvn_tlin  (:,:,:)
1131      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
1132      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
1133
1134      CALL div_cur_tan( nit000 )
1135
1136      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
1137      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
1138      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
1139      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
1140
1141      !--------------------------------------------------------------------
1142      ! Initialize the adjoint variables: dy^* = W dy
1143      !--------------------------------------------------------------------
1144      DO jk = 1, jpk
1145        DO jj = nldj, nlej
1146           DO ji = nldi, nlei
1147              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
1148                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1149                 &               * tmask(ji,jj,jk)
1150              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
1151                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1152                 &               * tmask(ji,jj,jk)
1153            END DO
1154         END DO
1155      END DO
1156      DO jk = 1, jpk
1157        DO jj = nldj, nlej
1158           DO ji = nldi, nlei
1159              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
1160                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1161              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
1162                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
1163            END DO
1164         END DO
1165      END DO
1166
1167      !--------------------------------------------------------------------
1168      ! Compute the scalar product: ( L dx )^T W dy
1169      !--------------------------------------------------------------------
1170
1171      zsp1_1 = DOT_PRODUCT( zhdivb_tlout, zhdivb_adin )
1172      zsp1_2 = DOT_PRODUCT( zhdivn_tlout, zhdivn_adin )
1173      zsp1_3 = DOT_PRODUCT( zrotb_tlout, zrotb_adin )
1174      zsp1_4 = DOT_PRODUCT( zrotn_tlout, zrotn_adin )
1175      zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4
1176
1177      !--------------------------------------------------------------------
1178      ! Call the adjoint routine: dx^* = L^T dy^*
1179      !--------------------------------------------------------------------
1180
1181      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
1182      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
1183      rotb_ad (:,:,:) = zrotb_adin(:,:,:)
1184      rotn_ad (:,:,:) = zrotn_adin(:,:,:)
1185
1186      CALL div_cur_adj( nit000 )
1187
1188      zun_adout  (:,:,:) = un_ad  (:,:,:)
1189      zvn_adout  (:,:,:) = vn_ad  (:,:,:)
1190      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
1191      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
1192
1193      !--------------------------------------------------------------------
1194      ! Compute the scalar product: dx^T L^T W dy
1195      !--------------------------------------------------------------------
1196
1197      zsp2_1 = DOT_PRODUCT( zun_tlin,   zun_adout   )
1198      zsp2_2 = DOT_PRODUCT( zvn_tlin,   zvn_adout   )
1199      zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
1200      zsp2_4 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
1201      zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
1202
1203      cl_name = 'div_cur_adj T3'
1204      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1205
1206
1207      DEALLOCATE( &
1208         & zun_tlin,     &
1209         & zvn_tlin,     &
1210         & zhdivn_tlin,  &
1211         & zrotn_tlin,   &
1212         & zhdivb_tlout, &
1213         & zhdivn_tlout, &
1214         & zrotb_tlout,  &
1215         & zrotn_tlout,  &
1216         & zhdivb_adin,  &
1217         & zhdivn_adin,  &
1218         & zrotb_adin,   &
1219         & zrotn_adin,   &
1220         & zun_adout,    &
1221         & zvn_adout,    &
1222         & zhdivn_adout, &
1223         & zrotn_adout,  &
1224         & znu,          &
1225         & znv           &
1226         & )
1227
1228   END SUBROUTINE div_cur_adj_tst
1229#endif
1230
1231   !!======================================================================
1232
1233END MODULE divcur_tam
Note: See TracBrowser for help on using the repository browser.