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

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/DYN/sshwzv_tam.F90 @ 5196

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

Correct a wrong adjoint computation in ssh_wzv_adj, see Ticket #1276
  • Property svn:executable set to *
File size: 36.1 KB
Line 
1MODULE sshwzv_tam
2#if defined key_tam
3   !!==============================================================================
4   !!                       ***  MODULE  sshwzv  ***
5   !! Ocean dynamics : sea surface height and vertical velocity
6   !!==============================================================================
7   !! History of the direct module:
8   !!            3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
9   !! History of the TAM module:
10   !!            3.2  !  2010-04  (F. Vigilant) Original code
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ssh_wzv        : after ssh & now vertical velocity
15   !!   ssh_nxt        : filter ans swap the ssh arrays
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE par_oce
19   USE in_out_manager
20   USE dom_oce
21   USE prtctl
22   USE phycst
23   USE lbclnk
24   USE lbclnk_tam
25   USE divcur_tam
26   USE cla_tam
27   USE oce_tam
28   USE sbc_oce_tam
29   USE gridrandom
30   USE dotprodfld
31   USE paresp
32   USE tstool_tam
33   USE lib_mpp
34   USE wrk_nemo
35   USE timing
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   ssh_wzv_tan       ! called by step.F90
41   PUBLIC   ssh_nxt_tan       ! called by step.F90
42   PUBLIC   ssh_wzv_adj       ! called by step.F90
43   PUBLIC   ssh_nxt_adj       ! called by step.F90
44   PUBLIC   ssh_wzv_adj_tst   ! called by tamtst.F90
45   PUBLIC   ssh_nxt_adj_tst   ! called by tamtst.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50
51CONTAINS
52
53   SUBROUTINE ssh_wzv_tan( kt , kdum )
54      !!----------------------------------------------------------------------
55      !!                ***  ROUTINE ssh_wzv_tan  ***
56      !!
57      !! ** Purpose of direct routine :
58      !!              compute the after ssh (ssha), the now vertical velocity
59      !!              and update the now vertical coordinate (lk_vvl=T).
60      !!
61      !! ** Method  : -
62      !!
63      !!              - Using the incompressibility hypothesis, the vertical
64      !!      velocity is computed by integrating the horizontal divergence
65      !!      from the bottom to the surface minus the scale factor evolution.
66      !!        The boundary conditions are w=0 at the bottom (no flux) and.
67      !!
68      !! ** action  :   ssha    : after sea surface height
69      !!                wn      : now vertical velocity
70      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
71      !!                          at u-, v-, f-point s
72      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
73      !!----------------------------------------------------------------------
74      !!
75      INTEGER, INTENT(in) ::   kt   ! time step
76      !!
77      INTEGER  ::   jk              ! dummy loop indices
78      REAL(wp) ::   z2dt, z1_rau0     ! temporary scalars
79      REAL(wp), POINTER, DIMENSION(:,:) ::   z2d, zhdivtl     ! 2D workspace
80      INTEGER, OPTIONAL            ::   kdum        ! dummy argument to compute only vertical velocity
81      !!----------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv_tan')
84      !
85      CALL wrk_alloc( jpi, jpj, z2d, zhdivtl )
86      z2d = 0._wp
87      zhdivtl = 0._wp
88      !
89      IF( kt == nit000 ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'ssh_wzv_tan : after sea surface height and now vertical velocity '
92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
93         !
94         wn_tl(:,:,jpk) = 0.0_wp                   ! bottom boundary condition: w=0 (set once for all)
95         !
96         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
97            IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
98            CALL abort
99         ENDIF
100         !
101      ENDIF
102      !                                           !------------------------------!
103      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
104                                                  !------------------------------!
105         IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
106         CALL abort
107         !
108      ENDIF
109
110      CALL div_cur_tan( kt )            ! Horizontal divergence & Relative vorticity
111
112      ! set time step size (Euler/Leapfrog)
113      z2dt = 2. * rdt
114      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
115
116      z1_rau0 = 0.5_wp / rau0
117
118      IF ( .NOT. PRESENT(kdum) ) THEN             ! jump ssh computing
119         !                                        !------------------------------!
120         !                                        !   After Sea Surface Height   !
121         !                                        !------------------------------!
122         zhdivtl(:,:) = 0.0_wp
123         DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
124            zhdivtl(:,:) = zhdivtl(:,:) + fse3t(:,:,jk) * hdivn_tl(:,:,jk)
125         END DO
126
127         !                                                ! Sea surface elevation time stepping
128         ssha_tl(:,:) = (  sshb_tl(:,:) - z2dt * ( z1_rau0 * ( emp_b_tl(:,:) + emp_tl(:,:) ) + zhdivtl(:,:) )  ) * tmask(:,:,1)
129
130         !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
131         IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
132            IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
133            CALL abort
134         ENDIF
135
136      ENDIF
137      !                                           !------------------------------!
138      !                                           !     Now Vertical Velocity    !
139      !                                           !------------------------------!
140      !                                                ! integrate from the bottom the hor. divergence
141      DO jk = jpkm1, 1, -1
142         wn_tl(:,:,jk) = wn_tl(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn_tl(:,:,jk)
143      END DO
144      !
145      CALL wrk_dealloc( jpi, jpj, z2d, zhdivtl )
146      !
147      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv_tan')
148      !
149   END SUBROUTINE ssh_wzv_tan
150
151   SUBROUTINE ssh_nxt_tan( kt )
152      !!----------------------------------------------------------------------
153      !!                    ***  ROUTINE ssh_nxt_tan  ***
154      !!
155      !! ** Purpose of the direct :
156      !!              achieve the sea surface  height time stepping by
157      !!              applying Asselin time filter and swapping the arrays
158      !!              ssha  already computed in ssh_wzv
159      !!
160      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
161      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
162      !!             sshn = ssha
163      !!
164      !! ** action  : - sshb, sshn   : before & now sea surface height
165      !!                               ready for the next time step
166      !!----------------------------------------------------------------------
167      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
168      !!
169      INTEGER  ::   ji, jj               ! dummy loop indices
170      REAL(wp) :: zec
171      !!----------------------------------------------------------------------
172      !
173      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt_tan')
174      !
175      IF( kt == nit000 ) THEN
176         IF(lwp) WRITE(numout,*)
177         IF(lwp) WRITE(numout,*) 'ssh_nxt_tan : next sea surface height (Asselin time filter + swap)'
178         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
179      ENDIF
180
181      ! Time filter and swap of the ssh
182      ! -------------------------------
183
184      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
185         !                   ! ---------------------- !
186         IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
187         CALL abort
188         !
189      ELSE                   ! fixed levels :   ssh at t-point only
190         !                   ! ------------ !
191         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
192            sshn_tl(:,:) = ssha_tl(:,:)                            ! now <-- after  (before already = now)
193            !
194         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
195            DO jj = 1, jpj
196               DO ji = 1, jpi                                ! before <-- now filtered
197                  sshb_tl(ji,jj) = sshn_tl(ji,jj) + atfp * ( sshb_tl(ji,jj) - 2 * sshn_tl(ji,jj) + ssha_tl(ji,jj) )
198                  sshn_tl(ji,jj) = ssha_tl(ji,jj)                  ! now <-- after
199               END DO
200            END DO
201         ENDIF
202      ENDIF
203      !
204#if defined key_agrif
205      ! Update velocity at AGRIF zoom boundaries
206      !IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn_tan( kt )
207      IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
208      CALL abort
209#endif
210      !
211      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt_tan')
212      !
213   END SUBROUTINE ssh_nxt_tan
214
215   SUBROUTINE ssh_wzv_adj( kt , kdum )
216      !!----------------------------------------------------------------------
217      !!                ***  ROUTINE ssh_wzv_adj  ***
218      !!
219      !! ** Purpose of direct routine :
220      !!              compute the after ssh (ssha), the now vertical velocity
221      !!              and update the now vertical coordinate (lk_vvl=T).
222      !!
223      !! ** Method  : -
224      !!
225      !!              - Using the incompressibility hypothesis, the vertical
226      !!      velocity is computed by integrating the horizontal divergence
227      !!      from the bottom to the surface minus the scale factor evolution.
228      !!        The boundary conditions are w=0 at the bottom (no flux) and.
229      !!
230      !! ** action  :   ssha    : after sea surface height
231      !!                wn      : now vertical velocity
232      !! if lk_vvl=T:   sshu_a, sshv_a, sshf_a  : after sea surface height
233      !!                          at u-, v-, f-point s
234      !!                hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
235      !!----------------------------------------------------------------------
236      !!
237      INTEGER, INTENT(in) ::   kt   ! time step
238      !!
239      INTEGER  ::   jk              ! dummy loop indices
240      REAL(wp) ::   z2dt, z1_rau0     ! temporary scalars
241      REAL(wp), POINTER, DIMENSION(:,:) ::  z2d, zhdivad     ! 2D workspace
242      INTEGER, OPTIONAL            ::   kdum        ! dummy argument to compute only vertical velocity
243      !!----------------------------------------------------------------------
244      !
245      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv_adj')
246      !
247      CALL wrk_alloc( jpi, jpj, z2d, zhdivad )
248      !
249      ! adjoint variable initialization
250      zhdivad = 0._wp
251      z2d = 0._wp
252
253      IF( kt == nitend ) THEN
254         IF(lwp) WRITE(numout,*)
255         IF(lwp) WRITE(numout,*) 'ssh_wzv_adj : after sea surface height and now vertical velocity '
256         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
257      ENDIF
258      !
259      !                                           !------------------------------!
260      IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
261                                                  !------------------------------!
262         IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
263         CALL abort
264         !
265      ENDIF
266      !                                           !------------------------------!
267      !                                           !     Now Vertical Velocity    !
268      !                                           !------------------------------!
269      !                                                ! integrate from the bottom the hor. divergence
270      DO jk = 1, jpkm1
271         hdivn_ad(:,:,jk  ) = hdivn_ad(:,:,jk  ) - wn_ad(:,:,jk) * fse3t_n(:,:,jk)
272         wn_ad(   :,:,jk+1) = wn_ad(   :,:,jk+1) + wn_ad(:,:,jk)
273         wn_ad(   :,:,jk  ) = 0.0_wp
274      END DO
275      !
276      ! set time step size (Euler/Leapfrog)
277      z2dt = 2. * rdt
278      IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
279
280      z1_rau0 = 0.5_wp / rau0
281
282      IF ( .NOT. PRESENT(kdum) ) THEN             ! jump ssh computing
283         !                                        !------------------------------!
284         !                                        !   After Sea Surface Height   !
285         !                                        !------------------------------!
286         !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
287         IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
288            IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
289            CALL abort
290         ENDIF
291         !                                                ! Sea surface elevation time stepping
292         ssha_ad( :,:) = ssha_ad( :,:) * tmask(:,:,1)
293         sshb_ad( :,:) = sshb_ad( :,:) + ssha_ad(:,:)
294         ssha_ad( :,:) = ssha_ad( :,:) * z2dt
295         zhdivad( :,:) = zhdivad( :,:) - ssha_ad(:,:)
296         ssha_ad( :,:) = ssha_ad( :,:) * z1_rau0 
297         emp_ad(  :,:) = emp_ad(  :,:) - ssha_ad(:,:)
298         emp_b_ad(:,:) = emp_b_ad(:,:) - ssha_ad(:,:)
299         ssha_ad(:,:) = 0.0_wp
300
301         DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
302            hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + fse3t(:,:,jk) * zhdivad(:,:)
303         END DO
304         zhdivad(:,:) = 0._wp
305
306      ENDIF
307
308      CALL div_cur_adj( kt )            ! Horizontal divergence & Relative vorticity
309      !
310      !
311      !
312      IF( kt == nit000 ) wn_ad(:,:,jpk) = 0._wp
313      !
314      !
315      !
316      CALL wrk_dealloc( jpi, jpj, z2d, zhdivad )
317      !
318      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv_adj')
319      !
320   END SUBROUTINE ssh_wzv_adj
321
322   SUBROUTINE ssh_nxt_adj( kt )
323      !!----------------------------------------------------------------------
324      !!                    ***  ROUTINE ssh_nxt_adj  ***
325      !!
326      !! ** Purpose of the direct :
327      !!              achieve the sea surface  height time stepping by
328      !!              applying Asselin time filter and swapping the arrays
329      !!              ssha  already computed in ssh_wzv
330      !!
331      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
332      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
333      !!             sshn = ssha
334      !!
335      !! ** action  : - sshb, sshn   : before & now sea surface height
336      !!                               ready for the next time step
337      !!----------------------------------------------------------------------
338      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
339      !!
340      INTEGER  ::   ji, jj               ! dummy loop indices
341      !!----------------------------------------------------------------------
342      !
343      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt_adj')
344      !
345
346      IF( kt == nitend ) THEN
347         IF(lwp) WRITE(numout,*)
348         IF(lwp) WRITE(numout,*) 'ssh_nxt_adj : next sea surface height (Asselin time filter + swap)'
349         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
350      ENDIF
351      ! Time filter and swap of the ssh
352      ! -------------------------------
353#if defined key_agrif
354      ! Update velocity at AGRIF zoom boundaries
355      !IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn_adj( kt )
356      IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
357      CALL abort
358#endif
359      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
360         !                   ! ---------------------- !
361         IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
362         CALL abort
363         !
364      ELSE                   ! fixed levels :   ssh at t-point only
365
366         !                   ! ------------ !
367         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
368            ssha_ad(:,:) = ssha_ad(:,:) + sshn_ad(:,:)
369            sshn_ad(:,:) = 0.0_wp
370            !
371         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
372            DO jj = jpj, 1, -1
373               DO ji = jpi, 1, -1                                ! before <-- now filtered
374                  ssha_ad(ji,jj) = ssha_ad(ji,jj) + sshn_ad(ji,jj)
375                  sshn_ad(ji,jj) = 0._wp
376                  sshn_ad(ji,jj) = sshn_ad(ji,jj) + (1.0_wp - 2.0 * atfp) * sshb_ad(ji,jj)
377                  ssha_ad(ji,jj) = ssha_ad(ji,jj) + atfp * sshb_ad(ji,jj)
378                  sshb_ad(ji,jj) = atfp * sshb_ad(ji,jj)
379               END DO
380            END DO
381         ENDIF
382      ENDIF
383      !
384      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt_adj')
385      !
386   END SUBROUTINE ssh_nxt_adj
387
388   SUBROUTINE ssh_wzv_adj_tst( kumadt )
389      !!-----------------------------------------------------------------------
390      !!
391      !!          ***  ROUTINE ssh_wzv_adj_tst : TEST OF wzv_adj  ***
392      !!
393      !! ** Purpose : Test the adjoint routine.
394      !!
395      !! ** Method  : Verify the scalar product
396      !!
397      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
398      !!
399      !!              where  L   = tangent routine
400      !!                     L^T = adjoint routine
401      !!                     W   = diagonal matrix of scale factors
402      !!                     dx  = input perturbation (random field)
403      !!                     dy  = L dx
404      !!
405      !! ** Action  : Separate tests are applied for the following dx and dy:
406      !!
407      !!            dx = ( un_tl, vn_tl, hdivn_tl, rotn_tl, emp_tl, sshb_tl ) and
408      !!            dy = ( hdivn_tl, hdivb_tl, rotn_tl, rotb_tl, wn_tl, ssha_tl )
409      !!
410      !! History :
411      !!        ! 2010-04 (F. Vigilant)
412      !!-----------------------------------------------------------------------
413
414      !! * Modules used
415      !! * Arguments
416      INTEGER, INTENT(IN) :: &
417         & kumadt             ! Output unit
418
419      !! * Local declarations
420      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
421         & zun_tlin,     & ! Tangent input: now u-velocity
422         & zvn_tlin,     & ! Tangent input: now v-velocity
423         & zhdivn_tlin,  & ! Tangent input: now horizontal divergence
424         & zrotn_tlin,   & ! Tangent input: now horizontal divergence
425         & zhdivn_tlout, & ! Tangent output: now horizontal divergence
426         & zrotn_tlout,  & ! Tangent output: now horizontal divergence
427         & zrotb_tlout,  & ! Tangent output: now horizontal divergence
428         & zhdivb_tlout, & ! Tangent output: now horizontal divergence
429         & zwn_tlout,    & ! Tangent output: now w-velocity
430         & zwn_adin,     & ! Adjoint input: now w-velocity
431         & zhdivn_adout, & ! Adjoint output: now horizontal divergence
432         & zrotn_adin,   & ! Adjoint input: now horizontal divergence
433         & zrotn_adout,  & ! Adjoint output: now horizontal divergence
434         & zrotb_adin,   & ! Adjoint input: now horizontal divergence
435         & zhdivn_adin,  & ! Adjoint input: now horizontal divergence
436         & zhdivb_adin,  & ! Adjoint output: now horizontal divergence
437         & zun_adout,    & ! Adjoint output: now horizontal divergence
438         & zvn_adout,    & ! Adjoint output: now horizontal divergence
439         & znu,          & ! 3D random field for u
440         & znv             ! 3D random field for v
441      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
442         & zsshb_tlin,   & ! Tangent input: before SSH
443         & zssha_tlout,  & ! Tangent input: before SSH
444         & zsshb_adout,  & ! Adjoint output: before SSH
445         & zssha_adin,   & ! Adjoint output: before SSH
446         & zemp_tlin,    & ! Tangent input: EmP
447         & zemp_adout,   & ! Adjoint output: EmP
448         & znssh,        & ! 2D random field for SSH
449         & znemp           ! 2D random field for EmP
450
451      INTEGER :: &
452         & ji,    &        ! dummy loop indices
453         & jj,    &
454         & jk
455      INTEGER, DIMENSION(jpi,jpj) :: &
456         & iseed_2d           ! 2D seed for the random number generator
457      REAL(KIND=wp) :: &
458                              ! random field standard deviation for:
459         & zstdssh,         & !   SSH
460         & zstdemp,         & !   EMP
461         & zsp1,            & ! scalar product involving the tangent routine
462         & zsp2,            & ! scalar product involving the adjoint routine
463         & zsp2_1,          & !   scalar product components
464         & zsp2_2,          &
465         & zsp2_3,          &
466         & zsp2_4,          &
467         & zsp2_5,          &
468         & zsp2_6,          &
469         & z2dt,            & ! temporary scalars
470         & zraur
471      CHARACTER (LEN=14) :: &
472         & cl_name
473
474      ! Allocate memory
475
476      ALLOCATE( &
477         & zhdivn_tlin(jpi,jpj,jpk),  &
478         & zhdivb_tlout(jpi,jpj,jpk), &
479         & zhdivn_tlout(jpi,jpj,jpk), &
480         & zrotn_tlin(jpi,jpj,jpk),   &
481         & zrotn_tlout(jpi,jpj,jpk),  &
482         & zrotb_tlout(jpi,jpj,jpk),  &
483         & zwn_tlout(jpi,jpj,jpk),    &
484         & zwn_adin(jpi,jpj,jpk),     &
485         & zhdivn_adout(jpi,jpj,jpk), &
486         & zhdivb_adin(jpi,jpj,jpk),  &
487         & zrotn_adin(jpi,jpj,jpk),   &
488         & zrotn_adout(jpi,jpj,jpk),  &
489         & zrotb_adin(jpi,jpj,jpk),   &
490         & zhdivn_adin(jpi,jpj,jpk),  &
491         & zun_tlin(jpi,jpj,jpk),     &
492         & zvn_tlin(jpi,jpj,jpk),     &
493         & zun_adout(jpi,jpj,jpk),    &
494         & zvn_adout(jpi,jpj,jpk),    &
495         & znu(jpi,jpj,jpk),          &
496         & znv(jpi,jpj,jpk)           &
497         & )
498      ALLOCATE( &
499         & zsshb_tlin(jpi,jpj),       &
500         & zsshb_adout(jpi,jpj),      &
501         & zssha_tlout(jpi,jpj),      &
502         & zssha_adin(jpi,jpj),       &
503         & zemp_tlin(jpi,jpj),        &
504         & zemp_adout(jpi,jpj),       &
505         & znssh(jpi,jpj),            &
506         & znemp(jpi,jpj)             &
507         & )
508
509
510      ! Initialize constants
511
512      z2dt  = 2.0_wp * rdt       ! time step: leap-frog
513      zraur = 1.0_wp / rau0      ! inverse density of pure water (m3/kg)
514
515      zhdivn_tlin(:,:,:) = 0.0_wp
516      zrotn_tlin(:,:,:) = 0.0_wp
517      zemp_tlin(:,:) = 0.0_wp
518      zsshb_tlin(:,:) = 0.0_wp
519      zun_tlin (:,:,:) = 0.0_wp
520      zvn_tlin (:,:,:) = 0.0_wp
521
522      zhdivn_tlout(:,:,:) = 0.0_wp
523      zhdivb_tlout(:,:,:) = 0.0_wp
524      zrotn_tlout(:,:,:)  = 0.0_wp
525      zrotb_tlout(:,:,:)  = 0.0_wp
526      zwn_tlout(:,:,:) = 0.0_wp
527      zssha_tlout(:,:) = 0.0_wp
528
529      zhdivn_adin(:,:,:) = 0.0_wp
530      zhdivb_adin(:,:,:) = 0.0_wp
531      zrotn_adin(:,:,:)  = 0.0_wp
532      zrotb_adin(:,:,:)  = 0.0_wp
533      zwn_adin(:,:,:) = 0.0_wp
534      zssha_adin(:,:) = 0.0_wp
535
536      zhdivn_adout(:,:,:) = 0.0_wp
537      zrotn_adout(:,:,:) = 0.0_wp
538      zemp_adout(:,:) = 0.0_wp
539      zsshb_adout(:,:) = 0.0_wp
540      zun_adout (:,:,:) = 0.0_wp
541      zvn_adout (:,:,:) = 0.0_wp
542
543      un_tl   (:,:,:) = 0.0_wp
544      vn_tl   (:,:,:) = 0.0_wp
545      hdivn_tl(:,:,:) = 0.0_wp
546      hdivb_tl(:,:,:) = 0.0_wp
547      rotn_tl (:,:,:) = 0.0_wp
548      rotb_tl (:,:,:) = 0.0_wp
549      wn_tl(:,:,:) = 0.0_wp
550      ssha_tl(:,:) = 0.0_wp
551      sshb_tl(:,:) = 0.0_wp
552      emp_tl(:,:) = 0.0_wp
553
554      un_ad   (:,:,:) = 0.0_wp
555      vn_ad   (:,:,:) = 0.0_wp
556      hdivn_ad(:,:,:) = 0.0_wp
557      hdivb_ad(:,:,:) = 0.0_wp
558      rotn_ad (:,:,:) = 0.0_wp
559      rotb_ad (:,:,:) = 0.0_wp
560      wn_ad(:,:,:) = 0.0_wp
561      sshb_ad(:,:) = 0.0_wp
562      ssha_ad(:,:) = 0.0_wp
563      emp_ad(:,:) = 0.0_wp
564
565      !=============================================================
566      ! 1) dx = ( un_tl, vn_tl, emp_tl, sshb_tl ) and dy = ( wn_tl )
567      !                   - or -
568      ! 2) dx = ( hdivn_tl ) and dy = ( wn_tl )
569      !=============================================================
570
571      !--------------------------------------------------------------------
572      ! Initialize the tangent input with random noise: dx
573      !--------------------------------------------------------------------
574
575      CALL grid_random(  znu, 'U', 0.0_wp, stdu )
576      CALL grid_random(  znv, 'V', 0.0_wp, stdv )
577
578      DO jk = 1, jpk
579         DO jj = nldj, nlej
580            DO ji = nldi, nlei
581               zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
582               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
583            END DO
584         END DO
585      END DO
586
587      CALL grid_random(  znssh, 'T', 0.0_wp, stdssh )
588      CALL grid_random(  znemp, 'T', 0.0_wp, stdssh )
589
590      DO jj = nldj, nlej
591         DO ji = nldi, nlei
592            zsshb_tlin(ji,jj) = znssh(ji,jj)
593            zemp_tlin (ji,jj) = znemp(ji,jj) / ( z2dt * zraur )
594         END DO
595      END DO
596
597      un_tl(:,:,:) = zun_tlin(:,:,:)
598      vn_tl(:,:,:) = zvn_tlin(:,:,:)
599      CALL div_cur_tan( nit000 )    ! Generate noise hdiv/rot fields
600
601      DO jk = 1, jpk
602         DO jj = nldj, nlej
603            DO ji = nldi, nlei
604               zhdivn_tlin(ji,jj,jk) = 0.5_wp * hdivn_tl(ji,jj,jk)
605               zrotn_tlin (ji,jj,jk) = 0.5_wp * rotn_tl (ji,jj,jk)
606            END DO
607         END DO
608      END DO
609
610      ! re-initialization to zero
611      un_tl   (:,:,:) = 0.0_wp
612      vn_tl   (:,:,:) = 0.0_wp
613      hdivb_tl(:,:,:) = 0.0_wp
614      hdivn_tl(:,:,:) = 0.0_wp
615      rotb_tl (:,:,:) = 0.0_wp
616      rotn_tl (:,:,:) = 0.0_wp
617
618      !--------------------------------------------------------------------
619      ! Call the tangent routine: dy = L dx
620      !--------------------------------------------------------------------
621
622      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
623      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
624      sshb_tl(:,:) = zsshb_tlin(:,:)
625      emp_tl (:,:) = zemp_tlin (:,:)
626      un_tl(:,:,:) = zun_tlin(:,:,:)
627      vn_tl(:,:,:) = zvn_tlin(:,:,:)
628
629      CALL ssh_wzv_tan( nit000+1 )
630
631      zwn_tlout(:,:,:) = wn_tl(:,:,:)
632      zssha_tlout(:,: ) = ssha_tl(:,:)
633      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
634      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
635      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
636      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
637      !--------------------------------------------------------------------
638      ! Initialize the adjoint variables: dy^* = W dy
639      !--------------------------------------------------------------------
640
641      DO jk = 1, jpk
642        DO jj = nldj, nlej
643           DO ji = nldi, nlei
644              zwn_adin(ji,jj,jk) = zwn_tlout(ji,jj,jk) &
645                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
646                 &               * tmask(ji,jj,jk)
647            END DO
648         END DO
649      END DO
650      DO jj = nldj, nlej
651         DO ji = nldi, nlei
652            zssha_adin(ji,jj) = zssha_tlout(ji,jj) &
653               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
654               &                   * tmask(ji,jj,1)
655         END DO
656      END DO
657      DO jk = 1, jpk
658        DO jj = nldj, nlej
659           DO ji = nldi, nlei
660              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
661                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
662                 &               * tmask(ji,jj,jk)
663              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
664                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
665                 &               * tmask(ji,jj,jk)
666            END DO
667         END DO
668      END DO
669      DO jk = 1, jpk
670        DO jj = nldj, nlej
671           DO ji = nldi, nlei
672              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
673                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
674              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
675                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
676            END DO
677         END DO
678      END DO
679
680      !--------------------------------------------------------------------
681      ! Compute the scalar product: ( L dx )^T W dy
682      !--------------------------------------------------------------------
683
684
685      zsp1 = DOT_PRODUCT( zwn_tlout, zwn_adin ) + DOT_PRODUCT( zssha_tlout, zssha_adin ) &
686           & +  DOT_PRODUCT( zhdivb_tlout, zhdivb_adin ) + DOT_PRODUCT( zhdivn_tlout, zhdivn_adin ) &
687           & +  DOT_PRODUCT( zrotb_tlout, zrotb_adin ) + DOT_PRODUCT( zrotn_tlout, zrotn_adin )
688      !--------------------------------------------------------------------
689      ! Call the adjoint routine: dx^* = L^T dy^*
690      !--------------------------------------------------------------------
691
692      wn_ad(:,:,:) = zwn_adin(:,:,:)
693      ssha_ad(:,:) = zssha_adin(:,:)
694      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
695      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
696      rotb_ad(:,:,:) = zrotb_adin(:,:,:)
697      rotn_ad(:,:,:) = zrotn_adin(:,:,:)
698
699      CALL ssh_wzv_adj( nit000+1 )
700
701      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
702      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
703      zsshb_adout(:,:) = sshb_ad(:,:)
704      zemp_adout (:,:) = emp_ad (:,:)
705      zun_adout(:,:,:) = un_ad(:,:,:)
706      zvn_adout(:,:,:) = vn_ad(:,:,:)
707
708      !--------------------------------------------------------------------
709      ! Compute the scalar product: dx^T L^T W dy
710      !--------------------------------------------------------------------
711
712      zsp2_1 = DOT_PRODUCT( zun_tlin,    zun_adout    )
713      zsp2_2 = DOT_PRODUCT( zvn_tlin,    zvn_adout    )
714      zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
715      zsp2_4 = DOT_PRODUCT( zemp_tlin,   zemp_adout   )
716      zsp2_5 = DOT_PRODUCT( zsshb_tlin,  zsshb_adout  )
717      zsp2_6 = DOT_PRODUCT( zrotn_tlin,  zrotn_adout  )
718
719      zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6
720
721      ! Compare the scalar products
722      ! 14 char:'12345678901234'
723      cl_name = 'sshwzv_adj    '
724      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
725
726   END SUBROUTINE ssh_wzv_adj_tst
727
728   SUBROUTINE ssh_nxt_adj_tst( kumadt )
729      !!-----------------------------------------------------------------------
730      !!
731      !!          ***  ROUTINE ssh_nxt_adj_tst : TEST OF nxt_adj  ***
732      !!
733      !! ** Purpose : Test the adjoint routine.
734      !!
735      !! ** Method  : Verify the scalar product
736      !!
737      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
738      !!
739      !!              where  L   = tangent routine
740      !!                     L^T = adjoint routine
741      !!                     W   = diagonal matrix of scale factors
742      !!                     dx  = input perturbation (random field)
743      !!                     dy  = L dx
744      !!
745      !! ** Action  : Separate tests are applied for the following dx and dy:
746      !!
747      !!            dx = ( sshb_tl, sshn_tl, ssha_tl ) and
748      !!            dy = ( ssb_tl, sshn_tl )
749      !!
750      !! History :
751      !!        ! 2010-05 (F. Vigilant)
752      !!-----------------------------------------------------------------------
753
754      !! * Modules used
755      !! * Arguments
756      INTEGER, INTENT(IN) :: &
757         & kumadt             ! Output unit
758
759      !! * Local declarations
760      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
761         & zsshb_tlin,   & ! Tangent input: before SSH
762         & zsshn_tlin,   & ! Tangent input: before SSH
763         & zssha_tlin,   & ! Tangent input: before SSH
764         & zsshb_tlout,  & ! Tangent output: before SSH
765         & zsshn_tlout,  & ! Tangent output: before SSH
766         & zsshb_adin,   & ! Adjoint input: before SSH
767         & zsshn_adin,   & ! Adjoint input: before SSH
768         & zsshb_adout,  & ! Adjoint output: before SSH
769         & zsshn_adout,  & ! Adjoint output: before SSH
770         & zssha_adout,  & ! Adjoint output: before SSH
771         & znssh           ! 2D random field for EmP
772
773      INTEGER :: &
774         & ji,    &        ! dummy loop indices
775         & jj,    &
776         & jk
777      INTEGER, DIMENSION(jpi,jpj) :: &
778         & iseed_2d           ! 2D seed for the random number generator
779      REAL(KIND=wp) :: &
780                              ! random field standard deviation for:
781         & zstdssh,         & !   SSH
782         & zsp1,            & ! scalar product involving the tangent routine
783         & zsp2,            & ! scalar product involving the adjoint routine
784         & zsp1_1,          & !   scalar product components
785         & zsp1_2,          &
786         & zsp2_1,          & !   scalar product components
787         & zsp2_2,          &
788         & zsp2_3,          &
789         & zsp2_4
790      CHARACTER (LEN=14) :: &
791         & cl_name
792
793      ! Allocate memory
794
795      ALLOCATE( &
796         & zsshb_tlin(jpi,jpj),       &
797         & zsshn_tlin(jpi,jpj),       &
798         & zssha_tlin(jpi,jpj),       &
799         & zsshb_tlout(jpi,jpj),      &
800         & zsshn_tlout(jpi,jpj),      &
801         & zsshb_adin(jpi,jpj),       &
802         & zsshn_adin(jpi,jpj),       &
803         & zsshb_adout(jpi,jpj),      &
804         & zsshn_adout(jpi,jpj),      &
805         & zssha_adout(jpi,jpj),      &
806         & znssh(jpi,jpj)             &
807         & )
808
809
810      ! Initialize constants
811
812      zsshb_tlin(:,:)  = 0.0_wp
813      zsshn_tlin(:,:)  = 0.0_wp
814      zssha_tlin(:,:)  = 0.0_wp
815
816      zsshb_tlout(:,:) = 0.0_wp
817      zsshn_tlout(:,:) = 0.0_wp
818
819      zsshb_adout(:,:) = 0.0_wp
820      zsshn_adout(:,:) = 0.0_wp
821      zssha_adout(:,:) = 0.0_wp
822
823      zsshb_adin(:,:)  = 0.0_wp
824      zsshn_adin(:,:)  = 0.0_wp
825
826      sshb_tl(:,:)     = 0.0_wp
827      sshn_tl(:,:)     = 0.0_wp
828      ssha_tl(:,:)     = 0.0_wp
829
830      sshb_ad(:,:)     = 0.0_wp
831      sshn_ad(:,:)     = 0.0_wp
832      ssha_ad(:,:)     = 0.0_wp
833
834      !=============================================================
835      ! dx = ( sshb_tl, sshn_tl, ssha_tl ) and dy = ( ssb_tl, sshn_tl )
836      !=============================================================
837
838      !--------------------------------------------------------------------
839      ! Initialize the tangent input with random noise: dx
840      !--------------------------------------------------------------------
841
842      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
843
844      DO jj = nldj, nlej
845         DO ji = nldi, nlei
846            zsshb_tlin(ji,jj) = znssh(ji,jj)
847         END DO
848      END DO
849
850      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
851
852      DO jj = nldj, nlej
853         DO ji = nldi, nlei
854            zsshn_tlin(ji,jj) = znssh(ji,jj)
855         END DO
856      END DO
857
858      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
859
860      DO jj = nldj, nlej
861         DO ji = nldi, nlei
862            zssha_tlin(ji,jj) = znssh(ji,jj)
863         END DO
864      END DO
865
866      !--------------------------------------------------------------------
867      ! Call the tangent routine: dy = L dx
868      !--------------------------------------------------------------------
869
870      sshb_tl(:,:) = zsshb_tlin(:,:)
871      sshn_tl(:,:) = zsshn_tlin(:,:)
872      ssha_tl(:,:) = zssha_tlin(:,:)
873
874      CALL ssh_nxt_tan( nit000+1 )
875
876      zsshb_tlout(:,: ) = sshb_tl(:,:)
877      zsshn_tlout(:,: ) = sshn_tl(:,:)
878      !--------------------------------------------------------------------
879      ! Initialize the adjoint variables: dy^* = W dy
880      !--------------------------------------------------------------------
881
882      DO jj = nldj, nlej
883         DO ji = nldi, nlei
884            zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) &
885               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
886               &                   * tmask(ji,jj,1)
887            zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
888               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
889               &                   * tmask(ji,jj,1)
890         END DO
891      END DO
892
893      !--------------------------------------------------------------------
894      ! Compute the scalar product: ( L dx )^T W dy
895      !--------------------------------------------------------------------
896      zsp1_1 = DOT_PRODUCT( zsshb_tlout, zsshb_adin )
897      zsp1_2 = DOT_PRODUCT( zsshn_tlout, zsshn_adin )
898
899      zsp1 = zsp1_1 + zsp1_2
900      !--------------------------------------------------------------------
901      ! Call the adjoint routine: dx^* = L^T dy^*
902      !--------------------------------------------------------------------
903
904      sshb_ad(:,:) = zsshb_adin(:,:)
905      sshn_ad(:,:) = zsshn_adin(:,:)
906
907      CALL ssh_nxt_adj( nit000+1 )
908
909      zsshb_adout(:,:) = sshb_ad(:,:)
910      zsshn_adout(:,:) = sshn_ad(:,:)
911      zssha_adout(:,:) = ssha_ad(:,:)
912
913      !--------------------------------------------------------------------
914      ! Compute the scalar product: dx^T L^T W dy
915      !--------------------------------------------------------------------
916
917      zsp2_1 = DOT_PRODUCT( zsshb_tlin,  zsshb_adout  )
918      zsp2_2 = DOT_PRODUCT( zsshn_tlin,  zsshn_adout  )
919      zsp2_3 = DOT_PRODUCT( zssha_tlin,  zssha_adout  )
920
921      zsp2 = zsp2_1 + zsp2_2 + zsp2_3
922
923      ! Compare the scalar products
924      ! 14 char:'12345678901234'
925      cl_name = 'sshnxt_adj    '
926      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
927
928   END SUBROUTINE ssh_nxt_adj_tst
929
930   !!======================================================================
931#endif
932
933END MODULE sshwzv_tam
Note: See TracBrowser for help on using the repository browser.