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

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/sshwzv_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: 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  ) - fse3t_n(:,:,jk) * wn_ad(:,:,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         sshb_ad(:,:) = sshb_ad(:,:) + ssha_ad(:,:)* tmask(:,:,1)
293         emp_ad( :,:) = emp_ad(:,:)  - z2dt * z1_rau0 * ssha_ad(:,:) * tmask(:,:,1)
294         emp_b_ad( :,:) = emp_b_ad(:,:) - z2dt * z1_rau0 * ssha_ad(:,:) * tmask(:,:,1)
295         zhdivad(:,:) = zhdivad(:,:) - z2dt * tmask(:,:,1) * ssha_ad(:,:)
296         ssha_ad(:,:) = 0.0_wp
297
298         DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
299            hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + fse3t(:,:,jk) * zhdivad(:,:)
300         END DO
301         zhdivad(:,:) = 0._wp
302
303      ENDIF
304
305      CALL div_cur_adj( kt )            ! Horizontal divergence & Relative vorticity
306      !
307      !
308      !
309      wn_ad(:,:,jpk) = 0._wp
310      !
311      !
312      !
313      CALL wrk_dealloc( jpi, jpj, z2d, zhdivad )
314      !
315      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv_adj')
316      !
317   END SUBROUTINE ssh_wzv_adj
318
319   SUBROUTINE ssh_nxt_adj( kt )
320      !!----------------------------------------------------------------------
321      !!                    ***  ROUTINE ssh_nxt_adj  ***
322      !!
323      !! ** Purpose of the direct :
324      !!              achieve the sea surface  height time stepping by
325      !!              applying Asselin time filter and swapping the arrays
326      !!              ssha  already computed in ssh_wzv
327      !!
328      !! ** Method  : - apply Asselin time fiter to now ssh and swap :
329      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha )
330      !!             sshn = ssha
331      !!
332      !! ** action  : - sshb, sshn   : before & now sea surface height
333      !!                               ready for the next time step
334      !!----------------------------------------------------------------------
335      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
336      !!
337      INTEGER  ::   ji, jj               ! dummy loop indices
338      !!----------------------------------------------------------------------
339      !
340      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt_adj')
341      !
342
343      IF( kt == nitend ) THEN
344         IF(lwp) WRITE(numout,*)
345         IF(lwp) WRITE(numout,*) 'ssh_nxt_adj : next sea surface height (Asselin time filter + swap)'
346         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
347      ENDIF
348      ! Time filter and swap of the ssh
349      ! -------------------------------
350#if defined key_agrif
351      ! Update velocity at AGRIF zoom boundaries
352      !IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn_adj( kt )
353      IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
354      CALL abort
355#endif
356      IF( lk_vvl ) THEN      ! Variable volume levels :   ssh at t-, u-, v, f-points
357         !                   ! ---------------------- !
358         IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
359         CALL abort
360         !
361      ELSE                   ! fixed levels :   ssh at t-point only
362
363         !                   ! ------------ !
364         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
365            ssha_ad(:,:) = ssha_ad(:,:) + sshn_ad(:,:)
366            sshn_ad(:,:) = 0.0_wp
367            !
368         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
369            DO jj = jpj, 1, -1
370               DO ji = jpi, 1, -1                                ! before <-- now filtered
371                  ssha_ad(ji,jj) = ssha_ad(ji,jj) + sshn_ad(ji,jj)
372                  sshn_ad(ji,jj) = 0._wp
373                  sshn_ad(ji,jj) = sshn_ad(ji,jj) + (1.0_wp - 2.0 * atfp) * sshb_ad(ji,jj)
374                  ssha_ad(ji,jj) = ssha_ad(ji,jj) + atfp * sshb_ad(ji,jj)
375                  sshb_ad(ji,jj) = atfp * sshb_ad(ji,jj)
376               END DO
377            END DO
378         ENDIF
379      ENDIF
380      !
381      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt_adj')
382      !
383   END SUBROUTINE ssh_nxt_adj
384
385   SUBROUTINE ssh_wzv_adj_tst( kumadt )
386      !!-----------------------------------------------------------------------
387      !!
388      !!          ***  ROUTINE ssh_wzv_adj_tst : TEST OF wzv_adj  ***
389      !!
390      !! ** Purpose : Test the adjoint routine.
391      !!
392      !! ** Method  : Verify the scalar product
393      !!
394      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
395      !!
396      !!              where  L   = tangent routine
397      !!                     L^T = adjoint routine
398      !!                     W   = diagonal matrix of scale factors
399      !!                     dx  = input perturbation (random field)
400      !!                     dy  = L dx
401      !!
402      !! ** Action  : Separate tests are applied for the following dx and dy:
403      !!
404      !!            dx = ( un_tl, vn_tl, hdivn_tl, rotn_tl, emp_tl, sshb_tl ) and
405      !!            dy = ( hdivn_tl, hdivb_tl, rotn_tl, rotb_tl, wn_tl, ssha_tl )
406      !!
407      !! History :
408      !!        ! 2010-04 (F. Vigilant)
409      !!-----------------------------------------------------------------------
410
411      !! * Modules used
412      !! * Arguments
413      INTEGER, INTENT(IN) :: &
414         & kumadt             ! Output unit
415
416      !! * Local declarations
417      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
418         & zun_tlin,     & ! Tangent input: now u-velocity
419         & zvn_tlin,     & ! Tangent input: now v-velocity
420         & zhdivn_tlin,  & ! Tangent input: now horizontal divergence
421         & zrotn_tlin,   & ! Tangent input: now horizontal divergence
422         & zhdivn_tlout, & ! Tangent output: now horizontal divergence
423         & zrotn_tlout,  & ! Tangent output: now horizontal divergence
424         & zrotb_tlout,  & ! Tangent output: now horizontal divergence
425         & zhdivb_tlout, & ! Tangent output: now horizontal divergence
426         & zwn_tlout,    & ! Tangent output: now w-velocity
427         & zwn_adin,     & ! Adjoint input: now w-velocity
428         & zhdivn_adout, & ! Adjoint output: now horizontal divergence
429         & zrotn_adin,   & ! Adjoint input: now horizontal divergence
430         & zrotn_adout,  & ! Adjoint output: now horizontal divergence
431         & zrotb_adin,   & ! Adjoint input: now horizontal divergence
432         & zhdivn_adin,  & ! Adjoint input: now horizontal divergence
433         & zhdivb_adin,  & ! Adjoint output: now horizontal divergence
434         & zun_adout,    & ! Adjoint output: now horizontal divergence
435         & zvn_adout,    & ! Adjoint output: now horizontal divergence
436         & znu,          & ! 3D random field for u
437         & znv             ! 3D random field for v
438      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
439         & zsshb_tlin,   & ! Tangent input: before SSH
440         & zssha_tlout,  & ! Tangent input: before SSH
441         & zsshb_adout,  & ! Adjoint output: before SSH
442         & zssha_adin,   & ! Adjoint output: before SSH
443         & zemp_tlin,    & ! Tangent input: EmP
444         & zemp_adout,   & ! Adjoint output: EmP
445         & znssh,        & ! 2D random field for SSH
446         & znemp           ! 2D random field for EmP
447
448      INTEGER :: &
449         & ji,    &        ! dummy loop indices
450         & jj,    &
451         & jk
452      INTEGER, DIMENSION(jpi,jpj) :: &
453         & iseed_2d           ! 2D seed for the random number generator
454      REAL(KIND=wp) :: &
455                              ! random field standard deviation for:
456         & zstdssh,         & !   SSH
457         & zstdemp,         & !   EMP
458         & zsp1,            & ! scalar product involving the tangent routine
459         & zsp2,            & ! scalar product involving the adjoint routine
460         & zsp2_1,          & !   scalar product components
461         & zsp2_2,          &
462         & zsp2_3,          &
463         & zsp2_4,          &
464         & zsp2_5,          &
465         & zsp2_6,          &
466         & z2dt,            & ! temporary scalars
467         & zraur
468      CHARACTER (LEN=14) :: &
469         & cl_name
470
471      ! Allocate memory
472
473      ALLOCATE( &
474         & zhdivn_tlin(jpi,jpj,jpk),  &
475         & zhdivb_tlout(jpi,jpj,jpk), &
476         & zhdivn_tlout(jpi,jpj,jpk), &
477         & zrotn_tlin(jpi,jpj,jpk),   &
478         & zrotn_tlout(jpi,jpj,jpk),  &
479         & zrotb_tlout(jpi,jpj,jpk),  &
480         & zwn_tlout(jpi,jpj,jpk),    &
481         & zwn_adin(jpi,jpj,jpk),     &
482         & zhdivn_adout(jpi,jpj,jpk), &
483         & zhdivb_adin(jpi,jpj,jpk),  &
484         & zrotn_adin(jpi,jpj,jpk),   &
485         & zrotn_adout(jpi,jpj,jpk),  &
486         & zrotb_adin(jpi,jpj,jpk),   &
487         & zhdivn_adin(jpi,jpj,jpk),  &
488         & zun_tlin(jpi,jpj,jpk),     &
489         & zvn_tlin(jpi,jpj,jpk),     &
490         & zun_adout(jpi,jpj,jpk),    &
491         & zvn_adout(jpi,jpj,jpk),    &
492         & znu(jpi,jpj,jpk),          &
493         & znv(jpi,jpj,jpk)           &
494         & )
495      ALLOCATE( &
496         & zsshb_tlin(jpi,jpj),       &
497         & zsshb_adout(jpi,jpj),      &
498         & zssha_tlout(jpi,jpj),      &
499         & zssha_adin(jpi,jpj),       &
500         & zemp_tlin(jpi,jpj),        &
501         & zemp_adout(jpi,jpj),       &
502         & znssh(jpi,jpj),            &
503         & znemp(jpi,jpj)             &
504         & )
505
506
507      ! Initialize constants
508
509      z2dt  = 2.0_wp * rdt       ! time step: leap-frog
510      zraur = 1.0_wp / rau0      ! inverse density of pure water (m3/kg)
511
512      zhdivn_tlin(:,:,:) = 0.0_wp
513      zrotn_tlin(:,:,:) = 0.0_wp
514      zemp_tlin(:,:) = 0.0_wp
515      zsshb_tlin(:,:) = 0.0_wp
516      zun_tlin (:,:,:) = 0.0_wp
517      zvn_tlin (:,:,:) = 0.0_wp
518
519      zhdivn_tlout(:,:,:) = 0.0_wp
520      zhdivb_tlout(:,:,:) = 0.0_wp
521      zrotn_tlout(:,:,:)  = 0.0_wp
522      zrotb_tlout(:,:,:)  = 0.0_wp
523      zwn_tlout(:,:,:) = 0.0_wp
524      zssha_tlout(:,:) = 0.0_wp
525
526      zhdivn_adin(:,:,:) = 0.0_wp
527      zhdivb_adin(:,:,:) = 0.0_wp
528      zrotn_adin(:,:,:)  = 0.0_wp
529      zrotb_adin(:,:,:)  = 0.0_wp
530      zwn_adin(:,:,:) = 0.0_wp
531      zssha_adin(:,:) = 0.0_wp
532
533      zhdivn_adout(:,:,:) = 0.0_wp
534      zrotn_adout(:,:,:) = 0.0_wp
535      zemp_adout(:,:) = 0.0_wp
536      zsshb_adout(:,:) = 0.0_wp
537      zun_adout (:,:,:) = 0.0_wp
538      zvn_adout (:,:,:) = 0.0_wp
539
540      un_tl   (:,:,:) = 0.0_wp
541      vn_tl   (:,:,:) = 0.0_wp
542      hdivn_tl(:,:,:) = 0.0_wp
543      hdivb_tl(:,:,:) = 0.0_wp
544      rotn_tl (:,:,:) = 0.0_wp
545      rotb_tl (:,:,:) = 0.0_wp
546      wn_tl(:,:,:) = 0.0_wp
547      ssha_tl(:,:) = 0.0_wp
548      sshb_tl(:,:) = 0.0_wp
549      emp_tl(:,:) = 0.0_wp
550
551      un_ad   (:,:,:) = 0.0_wp
552      vn_ad   (:,:,:) = 0.0_wp
553      hdivn_ad(:,:,:) = 0.0_wp
554      hdivb_ad(:,:,:) = 0.0_wp
555      rotn_ad (:,:,:) = 0.0_wp
556      rotb_ad (:,:,:) = 0.0_wp
557      wn_ad(:,:,:) = 0.0_wp
558      sshb_ad(:,:) = 0.0_wp
559      ssha_ad(:,:) = 0.0_wp
560      emp_ad(:,:) = 0.0_wp
561
562      !=============================================================
563      ! 1) dx = ( un_tl, vn_tl, emp_tl, sshb_tl ) and dy = ( wn_tl )
564      !                   - or -
565      ! 2) dx = ( hdivn_tl ) and dy = ( wn_tl )
566      !=============================================================
567
568      !--------------------------------------------------------------------
569      ! Initialize the tangent input with random noise: dx
570      !--------------------------------------------------------------------
571
572      CALL grid_random(  znu, 'U', 0.0_wp, stdu )
573      CALL grid_random(  znv, 'V', 0.0_wp, stdv )
574
575      DO jk = 1, jpk
576         DO jj = nldj, nlej
577            DO ji = nldi, nlei
578               zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
579               zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
580            END DO
581         END DO
582      END DO
583
584      CALL grid_random(  znssh, 'T', 0.0_wp, stdssh )
585      CALL grid_random(  znemp, 'T', 0.0_wp, stdssh )
586
587      DO jj = nldj, nlej
588         DO ji = nldi, nlei
589            zsshb_tlin(ji,jj) = znssh(ji,jj)
590            zemp_tlin (ji,jj) = znemp(ji,jj) / ( z2dt * zraur )
591         END DO
592      END DO
593
594      un_tl(:,:,:) = zun_tlin(:,:,:)
595      vn_tl(:,:,:) = zvn_tlin(:,:,:)
596      CALL div_cur_tan( nit000 )    ! Generate noise hdiv/rot fields
597
598      DO jk = 1, jpk
599         DO jj = nldj, nlej
600            DO ji = nldi, nlei
601               zhdivn_tlin(ji,jj,jk) = 0.5_wp * hdivn_tl(ji,jj,jk)
602               zrotn_tlin (ji,jj,jk) = 0.5_wp * rotn_tl (ji,jj,jk)
603            END DO
604         END DO
605      END DO
606
607      ! re-initialization to zero
608      un_tl   (:,:,:) = 0.0_wp
609      vn_tl   (:,:,:) = 0.0_wp
610      hdivb_tl(:,:,:) = 0.0_wp
611      hdivn_tl(:,:,:) = 0.0_wp
612      rotb_tl (:,:,:) = 0.0_wp
613      rotn_tl (:,:,:) = 0.0_wp
614
615      !--------------------------------------------------------------------
616      ! Call the tangent routine: dy = L dx
617      !--------------------------------------------------------------------
618
619      hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
620      rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
621      sshb_tl(:,:) = zsshb_tlin(:,:)
622      emp_tl (:,:) = zemp_tlin (:,:)
623      un_tl(:,:,:) = zun_tlin(:,:,:)
624      vn_tl(:,:,:) = zvn_tlin(:,:,:)
625
626      CALL ssh_wzv_tan( nit000+1 )
627
628      zwn_tlout(:,:,:) = wn_tl(:,:,:)
629      zssha_tlout(:,: ) = ssha_tl(:,:)
630      zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
631      zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
632      zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
633      zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
634      !--------------------------------------------------------------------
635      ! Initialize the adjoint variables: dy^* = W dy
636      !--------------------------------------------------------------------
637
638      DO jk = 1, jpk
639        DO jj = nldj, nlej
640           DO ji = nldi, nlei
641              zwn_adin(ji,jj,jk) = zwn_tlout(ji,jj,jk) &
642                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
643                 &               * tmask(ji,jj,jk)
644            END DO
645         END DO
646      END DO
647      DO jj = nldj, nlej
648         DO ji = nldi, nlei
649            zssha_adin(ji,jj) = zssha_tlout(ji,jj) &
650               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
651               &                   * tmask(ji,jj,1)
652         END DO
653      END DO
654      DO jk = 1, jpk
655        DO jj = nldj, nlej
656           DO ji = nldi, nlei
657              zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
658                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
659                 &               * tmask(ji,jj,jk)
660              zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
661                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
662                 &               * tmask(ji,jj,jk)
663            END DO
664         END DO
665      END DO
666      DO jk = 1, jpk
667        DO jj = nldj, nlej
668           DO ji = nldi, nlei
669              zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
670                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
671              zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
672                 &               * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
673            END DO
674         END DO
675      END DO
676
677      !--------------------------------------------------------------------
678      ! Compute the scalar product: ( L dx )^T W dy
679      !--------------------------------------------------------------------
680
681
682      zsp1 = DOT_PRODUCT( zwn_tlout, zwn_adin ) + DOT_PRODUCT( zssha_tlout, zssha_adin ) &
683           & +  DOT_PRODUCT( zhdivb_tlout, zhdivb_adin ) + DOT_PRODUCT( zhdivn_tlout, zhdivn_adin ) &
684           & +  DOT_PRODUCT( zrotb_tlout, zrotb_adin ) + DOT_PRODUCT( zrotn_tlout, zrotn_adin )
685      !--------------------------------------------------------------------
686      ! Call the adjoint routine: dx^* = L^T dy^*
687      !--------------------------------------------------------------------
688
689      wn_ad(:,:,:) = zwn_adin(:,:,:)
690      ssha_ad(:,:) = zssha_adin(:,:)
691      hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
692      hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
693      rotb_ad(:,:,:) = zrotb_adin(:,:,:)
694      rotn_ad(:,:,:) = zrotn_adin(:,:,:)
695
696      CALL ssh_wzv_adj( nit000+1 )
697
698      zrotn_adout(:,:,:) = rotn_ad(:,:,:)
699      zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
700      zsshb_adout(:,:) = sshb_ad(:,:)
701      zemp_adout (:,:) = emp_ad (:,:)
702      zun_adout(:,:,:) = un_ad(:,:,:)
703      zvn_adout(:,:,:) = vn_ad(:,:,:)
704
705      !--------------------------------------------------------------------
706      ! Compute the scalar product: dx^T L^T W dy
707      !--------------------------------------------------------------------
708
709      zsp2_1 = DOT_PRODUCT( zun_tlin,    zun_adout    )
710      zsp2_2 = DOT_PRODUCT( zvn_tlin,    zvn_adout    )
711      zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
712      zsp2_4 = DOT_PRODUCT( zemp_tlin,   zemp_adout   )
713      zsp2_5 = DOT_PRODUCT( zsshb_tlin,  zsshb_adout  )
714      zsp2_6 = DOT_PRODUCT( zrotn_tlin,  zrotn_adout  )
715
716      zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6
717
718      ! Compare the scalar products
719      ! 14 char:'12345678901234'
720      cl_name = 'sshwzv_adj    '
721      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
722
723   END SUBROUTINE ssh_wzv_adj_tst
724
725   SUBROUTINE ssh_nxt_adj_tst( kumadt )
726      !!-----------------------------------------------------------------------
727      !!
728      !!          ***  ROUTINE ssh_nxt_adj_tst : TEST OF nxt_adj  ***
729      !!
730      !! ** Purpose : Test the adjoint routine.
731      !!
732      !! ** Method  : Verify the scalar product
733      !!
734      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
735      !!
736      !!              where  L   = tangent routine
737      !!                     L^T = adjoint routine
738      !!                     W   = diagonal matrix of scale factors
739      !!                     dx  = input perturbation (random field)
740      !!                     dy  = L dx
741      !!
742      !! ** Action  : Separate tests are applied for the following dx and dy:
743      !!
744      !!            dx = ( sshb_tl, sshn_tl, ssha_tl ) and
745      !!            dy = ( ssb_tl, sshn_tl )
746      !!
747      !! History :
748      !!        ! 2010-05 (F. Vigilant)
749      !!-----------------------------------------------------------------------
750
751      !! * Modules used
752      !! * Arguments
753      INTEGER, INTENT(IN) :: &
754         & kumadt             ! Output unit
755
756      !! * Local declarations
757      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
758         & zsshb_tlin,   & ! Tangent input: before SSH
759         & zsshn_tlin,   & ! Tangent input: before SSH
760         & zssha_tlin,   & ! Tangent input: before SSH
761         & zsshb_tlout,  & ! Tangent output: before SSH
762         & zsshn_tlout,  & ! Tangent output: before SSH
763         & zsshb_adin,   & ! Adjoint input: before SSH
764         & zsshn_adin,   & ! Adjoint input: before SSH
765         & zsshb_adout,  & ! Adjoint output: before SSH
766         & zsshn_adout,  & ! Adjoint output: before SSH
767         & zssha_adout,  & ! Adjoint output: before SSH
768         & znssh           ! 2D random field for EmP
769
770      INTEGER :: &
771         & ji,    &        ! dummy loop indices
772         & jj,    &
773         & jk
774      INTEGER, DIMENSION(jpi,jpj) :: &
775         & iseed_2d           ! 2D seed for the random number generator
776      REAL(KIND=wp) :: &
777                              ! random field standard deviation for:
778         & zstdssh,         & !   SSH
779         & zsp1,            & ! scalar product involving the tangent routine
780         & zsp2,            & ! scalar product involving the adjoint routine
781         & zsp1_1,          & !   scalar product components
782         & zsp1_2,          &
783         & zsp2_1,          & !   scalar product components
784         & zsp2_2,          &
785         & zsp2_3,          &
786         & zsp2_4
787      CHARACTER (LEN=14) :: &
788         & cl_name
789
790      ! Allocate memory
791
792      ALLOCATE( &
793         & zsshb_tlin(jpi,jpj),       &
794         & zsshn_tlin(jpi,jpj),       &
795         & zssha_tlin(jpi,jpj),       &
796         & zsshb_tlout(jpi,jpj),      &
797         & zsshn_tlout(jpi,jpj),      &
798         & zsshb_adin(jpi,jpj),       &
799         & zsshn_adin(jpi,jpj),       &
800         & zsshb_adout(jpi,jpj),      &
801         & zsshn_adout(jpi,jpj),      &
802         & zssha_adout(jpi,jpj),      &
803         & znssh(jpi,jpj)             &
804         & )
805
806
807      ! Initialize constants
808
809      zsshb_tlin(:,:)  = 0.0_wp
810      zsshn_tlin(:,:)  = 0.0_wp
811      zssha_tlin(:,:)  = 0.0_wp
812
813      zsshb_tlout(:,:) = 0.0_wp
814      zsshn_tlout(:,:) = 0.0_wp
815
816      zsshb_adout(:,:) = 0.0_wp
817      zsshn_adout(:,:) = 0.0_wp
818      zssha_adout(:,:) = 0.0_wp
819
820      zsshb_adin(:,:)  = 0.0_wp
821      zsshn_adin(:,:)  = 0.0_wp
822
823      sshb_tl(:,:)     = 0.0_wp
824      sshn_tl(:,:)     = 0.0_wp
825      ssha_tl(:,:)     = 0.0_wp
826
827      sshb_ad(:,:)     = 0.0_wp
828      sshn_ad(:,:)     = 0.0_wp
829      ssha_ad(:,:)     = 0.0_wp
830
831      !=============================================================
832      ! dx = ( sshb_tl, sshn_tl, ssha_tl ) and dy = ( ssb_tl, sshn_tl )
833      !=============================================================
834
835      !--------------------------------------------------------------------
836      ! Initialize the tangent input with random noise: dx
837      !--------------------------------------------------------------------
838
839      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
840
841      DO jj = nldj, nlej
842         DO ji = nldi, nlei
843            zsshb_tlin(ji,jj) = znssh(ji,jj)
844         END DO
845      END DO
846
847      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
848
849      DO jj = nldj, nlej
850         DO ji = nldi, nlei
851            zsshn_tlin(ji,jj) = znssh(ji,jj)
852         END DO
853      END DO
854
855      CALL grid_random( znssh, 'T', 0.0_wp, stdssh )
856
857      DO jj = nldj, nlej
858         DO ji = nldi, nlei
859            zssha_tlin(ji,jj) = znssh(ji,jj)
860         END DO
861      END DO
862
863      !--------------------------------------------------------------------
864      ! Call the tangent routine: dy = L dx
865      !--------------------------------------------------------------------
866
867      sshb_tl(:,:) = zsshb_tlin(:,:)
868      sshn_tl(:,:) = zsshn_tlin(:,:)
869      ssha_tl(:,:) = zssha_tlin(:,:)
870
871      CALL ssh_nxt_tan( nit000+1 )
872
873      zsshb_tlout(:,: ) = sshb_tl(:,:)
874      zsshn_tlout(:,: ) = sshn_tl(:,:)
875      !--------------------------------------------------------------------
876      ! Initialize the adjoint variables: dy^* = W dy
877      !--------------------------------------------------------------------
878
879      DO jj = nldj, nlej
880         DO ji = nldi, nlei
881            zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) &
882               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
883               &                   * tmask(ji,jj,1)
884            zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
885               &                   * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
886               &                   * tmask(ji,jj,1)
887         END DO
888      END DO
889
890      !--------------------------------------------------------------------
891      ! Compute the scalar product: ( L dx )^T W dy
892      !--------------------------------------------------------------------
893      zsp1_1 = DOT_PRODUCT( zsshb_tlout, zsshb_adin )
894      zsp1_2 = DOT_PRODUCT( zsshn_tlout, zsshn_adin )
895
896      zsp1 = zsp1_1 + zsp1_2
897      !--------------------------------------------------------------------
898      ! Call the adjoint routine: dx^* = L^T dy^*
899      !--------------------------------------------------------------------
900
901      sshb_ad(:,:) = zsshb_adin(:,:)
902      sshn_ad(:,:) = zsshn_adin(:,:)
903
904      CALL ssh_nxt_adj( nit000+1 )
905
906      zsshb_adout(:,:) = sshb_ad(:,:)
907      zsshn_adout(:,:) = sshn_ad(:,:)
908      zssha_adout(:,:) = ssha_ad(:,:)
909
910      !--------------------------------------------------------------------
911      ! Compute the scalar product: dx^T L^T W dy
912      !--------------------------------------------------------------------
913
914      zsp2_1 = DOT_PRODUCT( zsshb_tlin,  zsshb_adout  )
915      zsp2_2 = DOT_PRODUCT( zsshn_tlin,  zsshn_adout  )
916      zsp2_3 = DOT_PRODUCT( zssha_tlin,  zssha_adout  )
917
918      zsp2 = zsp2_1 + zsp2_2 + zsp2_3
919
920      ! Compare the scalar products
921      ! 14 char:'12345678901234'
922      cl_name = 'sshnxt_adj    '
923      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
924
925   END SUBROUTINE ssh_nxt_adj_tst
926
927   !!======================================================================
928#endif
929
930END MODULE sshwzv_tam
Note: See TracBrowser for help on using the repository browser.