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

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

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