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.
dynspg_flt_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/dynspg_flt_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: 41.5 KB
Line 
1MODULE dynspg_flt_tam
2   !!----------------------------------------------------------------------
3   !!    This software is governed by the CeCILL licence (Version 2)
4   !!----------------------------------------------------------------------
5#if defined key_tam
6   !!======================================================================
7   !!  ***  MODULE  dynspg_flt_tam : TANGENT/ADJOINT OF MODULE dynspg_flt  ***
8   !!
9   !! Ocean dynamics:  surface pressure gradient trend
10   !!
11   !!======================================================================
12   !! History of the direct module:
13   !! History    OPA  !  1998-05  (G. Roullet)  free surface
14   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2
15   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module
16   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
17   !!            1.0  !  2004-08  (C. Talandier) New trends organization
18   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
19   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
20   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.
21   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
22   !! History of the TAM module:
23   !!            9.0  ! 2008-08  (A. Vidard) skeleton
24   !!             -   ! 2008-08  (A. Weaver) original version
25   !!   NEMO     3.2  ! 2010-04  (F. Vigilant) converison to 3.2
26   !!   NEMO     3.4  ! 2012-07  (P.-A. Bouttier) converison to 3.4
27   !!----------------------------------------------------------------------
28# if defined key_dynspg_flt   ||   defined key_esopa
29   !!----------------------------------------------------------------------
30   !!   'key_dynspg_flt'                              filtered free surface
31   !!----------------------------------------------------------------------
32   !!----------------------------------------------------------------------
33   !!   dyn_spg_flt_tan  : update the momentum trend with the surface pressure
34   !!                      gradient in the filtered free surface case
35   !!                      (tangent routine)
36   !!   dyn_spg_flt_adj  : update the momentum trend with the surface pressure
37   !!                      gradient in the filtered free surface case
38   !!                      (adjoint routine)
39   !!   dyn_spg_flt_adj_tst : Test of the adjoint routine
40   !!----------------------------------------------------------------------
41   USE prtctl                  ! Print control
42   USE par_oce
43   USE in_out_manager
44   USE phycst
45   USE lib_mpp
46   USE dom_oce
47   USE solver
48   USE sol_oce
49   USE oce_tam
50   USE sbc_oce_tam
51   USE sol_oce
52   USE sol_oce_tam
53   USE solsor_tam
54   USE lbclnk
55   USE lbclnk_tam
56   USE gridrandom
57   USE dotprodfld
58   USE paresp
59   USE dynadv
60   USE cla_tam
61   USE tstool_tam
62   USE wrk_nemo        ! Memory Allocation
63   USE lib_fortran
64   USE timing
65
66
67
68   IMPLICIT NONE
69   PRIVATE
70
71   !! * Accessibility
72   PUBLIC dyn_spg_flt_tan,     & ! routine called by step_tan.F90
73      &   dyn_spg_flt_adj,     & ! routine called by step_adj.F90
74      &   dyn_spg_flt_adj_tst    ! routine called by the tst.F90
75   !! * Substitutions
76#  include "domzgr_substitute.h90"
77#  include "vectopt_loop_substitute.h90"
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   SUBROUTINE dyn_spg_flt_tan( kt, kindic )
82      !!---------------------------------------------------------------------
83      !!                  ***  routine dyn_spg_flt_tan  ***
84      !!
85      !! ** Purpose of the direct routine:
86      !!      Compute the now trend due to the surface pressure
87      !!      gradient in case of filtered free surface formulation  and add
88      !!      it to the general trend of momentum equation.
89      !!
90      !! ** Method of the direct routine:
91      !!      Filtered free surface formulation. The surface
92      !!      pressure gradient is given by:
93      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
94      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
95      !!      where sshn is the free surface elevation and btda is the after
96      !!      time derivative of the free surface elevation
97      !!       -1- evaluate the surface presure trend (including the addi-
98      !!      tional force) in three steps:
99      !!        a- compute the right hand side of the elliptic equation:
100      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
101      !!         where (spgu,spgv) are given by:
102      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
103      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
104      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
105      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
106      !!         and define the first guess from previous computation :
107      !!            zbtd = btda
108      !!            btda = 2 zbtd - btdb
109      !!            btdb = zbtd
110      !!        b- compute the relative accuracy to be reached by the
111      !!         iterative solver
112      !!        c- apply the solver by a call to sol... routine
113      !!       -2- compute and add the free surface pressure gradient inclu-
114      !!      ding the additional force used to stabilize the equation.
115      !!
116      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
117      !!
118      !! References : Roullet and Madec 1999, JGR.
119      !!---------------------------------------------------------------------
120      !! * Arguments
121      INTEGER, INTENT( IN  ) :: &
122         &  kt         ! ocean time-step index
123      INTEGER, INTENT( OUT ) :: &
124         &  kindic     ! solver convergence flag (<0 if not converge)
125
126      !! * Local declarations
127      INTEGER  :: &
128         & ji,    &     ! dummy loop indices
129         & jj,    &
130         & jk
131      REAL(wp) ::  &
132         & z2dt,   & ! temporary scalars
133         & z2dtg,  &
134         & zgcb,   &
135         & zbtd,   &
136         & ztdgu,  &
137         & ztdgv
138      !!----------------------------------------------------------------------
139      !
140      !
141      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt_tan')
142      !
143      IF( kt == nit000 ) THEN
144
145         IF(lwp) WRITE(numout,*)
146         IF(lwp) WRITE(numout,*) 'dyn_spg_flt_tan : surface pressure gradient trend'
147         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   (free surface constant volume case)'
148         ! set to zero free surface specific arrays
149         spgu_tl(:,:) = 0.0_wp                    ! surface pressure gradient (i-direction)
150         spgv_tl(:,:) = 0.0_wp                    ! surface pressure gradient (j-direction)
151         !CALL solver_init( nit000 )               ! Elliptic solver initialisation
152      ENDIF
153      ! Local constant initialization
154      z2dt = 2. * rdt                                             ! time step: leap-frog
155      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
156      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
157      z2dtg  = grav * z2dt
158      ! Evaluate the masked next velocity (effect of the additional force not included)
159      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
160         CALL ctl_stop('key_vvl is not implemented in TAM yet')
161         !
162      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
163         !
164         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               spgu_tl(ji,jj) = - grav * ( sshn_tl(ji+1,jj) - sshn_tl(ji,jj) ) / e1u(ji,jj)
167               spgv_tl(ji,jj) = - grav * ( sshn_tl(ji,jj+1) - sshn_tl(ji,jj) ) / e2v(ji,jj)
168            END DO
169         END DO
170         DO jk = 1, jpkm1              ! unweighted time stepping
171            DO jj = 2, jpjm1
172               DO ji = fs_2, fs_jpim1   ! vector opt.
173                  ua_tl(ji,jj,jk) = (  ub_tl(ji,jj,jk) + z2dt * ( ua_tl(ji,jj,jk) + spgu_tl(ji,jj) )  ) * umask(ji,jj,jk)
174                  va_tl(ji,jj,jk) = (  vb_tl(ji,jj,jk) + z2dt * ( va_tl(ji,jj,jk) + spgv_tl(ji,jj) )  ) * vmask(ji,jj,jk)
175               END DO
176            END DO
177         END DO
178         !
179      ENDIF
180      IF( nn_cla == 1 )   CALL cla_dynspg_tan( kt )      ! Cross Land Advection (update (ua,va))
181
182      ! compute the next vertically averaged velocity (effect of the additional force not included)
183      ! ---------------------------------------------
184      DO jj = 2, jpjm1
185         DO ji = fs_2, fs_jpim1   ! vector opt.
186            spgu_tl(ji,jj) = 0.0_wp
187            spgv_tl(ji,jj) = 0.0_wp
188         END DO
189      END DO
190
191      ! vertical sum
192!CDIR NOLOOPCHG
193      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
194         DO jk = 1, jpkm1
195            DO ji = 1, jpij
196               spgu_tl(ji,1) = spgu_tl(ji,1) + fse3u(ji,1,jk) * ua_tl(ji,1,jk)
197               spgv_tl(ji,1) = spgv_tl(ji,1) + fse3v(ji,1,jk) * va_tl(ji,1,jk)
198            END DO
199         END DO
200      ELSE                        ! No  vector opt.
201         DO jk = 1, jpkm1
202            DO jj = 2, jpjm1
203               DO ji = 2, jpim1
204                  spgu_tl(ji,jj) = spgu_tl(ji,jj) + fse3u(ji,jj,jk) * ua_tl(ji,jj,jk)
205                  spgv_tl(ji,jj) = spgv_tl(ji,jj) + fse3v(ji,jj,jk) * va_tl(ji,jj,jk)
206               END DO
207            END DO
208         END DO
209      ENDIF
210
211      ! transport: multiplied by the horizontal scale factor
212      DO jj = 2, jpjm1
213         DO ji = fs_2, fs_jpim1   ! vector opt.
214            spgu_tl(ji,jj) = spgu_tl(ji,jj) * e2u(ji,jj)
215            spgv_tl(ji,jj) = spgv_tl(ji,jj) * e1v(ji,jj)
216         END DO
217      END DO
218
219      CALL lbc_lnk( spgu_tl, 'U', -1.0_wp )       ! lateral boundary conditions
220      CALL lbc_lnk( spgv_tl, 'V', -1.0_wp )
221
222      IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_tan: lk_vvl is not available' )
223
224      ! Right hand side of the elliptic equation and first guess
225      ! -----------------------------------------------------------
226      DO jj = 2, jpjm1
227         DO ji = fs_2, fs_jpim1   ! vector opt.
228            ! Divergence of the after vertically averaged velocity
229            zgcb =  spgu_tl(ji,jj) - spgu_tl(ji-1,jj)   &
230               &  + spgv_tl(ji,jj) - spgv_tl(ji,jj-1)
231            gcb_tl(ji,jj) = gcdprc(ji,jj) * zgcb
232            ! First guess of the after barotropic transport divergence
233            zbtd = gcx_tl(ji,jj)
234            gcx_tl (ji,jj) = 2.0_wp * zbtd - gcxb_tl(ji,jj)
235            gcxb_tl(ji,jj) =          zbtd
236         END DO
237      END DO
238      ! apply the lateral boundary conditions
239      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb_tl, c_solver_pt, 1.0_wp )
240
241      ! Relative precision
242      ! ------------------
243
244      rnorme = GLOB_SUM( gcb_tl(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb_tl(1:jpi,1:jpj) * bmask(:,:) )
245
246      epsr = eps * eps * rnorme
247      ncut = 0
248      ! if rnorme is 0, the solution is 0, the solver is not called
249      IF( rnorme == 0.0_wp ) THEN
250         gcx_tl(:,:) = 0.0_wp
251         res   = 0.0_wp
252         niter = 0
253         ncut  = 999
254      ENDIF
255
256      ! Evaluate the next transport divergence
257      ! --------------------------------------
258      !    Iterarive solver for the elliptic equation (except IF sol.=0)
259      !    (output in gcx with boundary conditions applied)
260      kindic = 0
261      IF( ncut == 0 ) THEN
262         IF    ( nn_solv == 1 ) THEN  ;  CALL ctl_stop('sol_pcg_tan not available in TAM yet')   ! diagonal preconditioned conjuguate gradient
263         ELSEIF( nn_solv == 2 ) THEN  ;  CALL sol_sor_tan( kt, kindic )   ! successive-over-relaxation
264         ENDIF
265      ENDIF
266
267      ! Transport divergence gradient multiplied by z2dt
268      ! --------------------------------------------====
269      DO jj = 2, jpjm1
270         DO ji = fs_2, fs_jpim1   ! vector opt.
271            ! trend of Transport divergence gradient
272            ztdgu = z2dtg * ( gcx_tl(ji+1,jj  ) - gcx_tl(ji,jj) ) / e1u(ji,jj)
273            ztdgv = z2dtg * ( gcx_tl(ji  ,jj+1) - gcx_tl(ji,jj) ) / e2v(ji,jj)
274            ! multiplied by z2dt
275            spgu_tl(ji,jj) = z2dt * ztdgu
276            spgv_tl(ji,jj) = z2dt * ztdgv
277         END DO
278      END DO
279
280      ! Add the trends multiplied by z2dt to the after velocity
281      ! -------------------------------------------------------
282      !     ( c a u t i o n : (ua_tl,va_tl) here are the after velocity not the
283      !                       trend, the leap-frog time stepping will not
284      !                       be done in dynnxt_tan.F90 routine)
285      DO jk = 1, jpkm1
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1   ! vector opt.
288               ua_tl(ji,jj,jk) = ( ua_tl(ji,jj,jk) + spgu_tl(ji,jj) ) * umask(ji,jj,jk)
289               va_tl(ji,jj,jk) = ( va_tl(ji,jj,jk) + spgv_tl(ji,jj) ) * vmask(ji,jj,jk)
290            END DO
291         END DO
292      END DO
293      !
294      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt_tan')
295      !
296   END SUBROUTINE dyn_spg_flt_tan
297
298   SUBROUTINE dyn_spg_flt_adj( kt, kindic )
299      !!----------------------------------------------------------------------
300      !!                  ***  routine dyn_spg_flt_adj  ***
301      !!
302      !! ** Purpose of the direct routine:
303      !!      Compute the now trend due to the surface pressure
304      !!      gradient in case of filtered free surface formulation  and add
305      !!      it to the general trend of momentum equation.
306      !!
307      !! ** Method of the direct routine:
308      !!      Filtered free surface formulation. The surface
309      !!      pressure gradient is given by:
310      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
311      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
312      !!      where sshn is the free surface elevation and btda is the after
313      !!      time derivative of the free surface elevation
314      !!       -1- evaluate the surface presure trend (including the addi-
315      !!      tional force) in three steps:
316      !!        a- compute the right hand side of the elliptic equation:
317      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
318      !!         where (spgu,spgv) are given by:
319      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
320      !!                 - grav 2 rdt hu /e1u di[sshn + emp]
321      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
322      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
323      !!         and define the first guess from previous computation :
324      !!            zbtd = btda
325      !!            btda = 2 zbtd - btdb
326      !!            btdb = zbtd
327      !!        b- compute the relative accuracy to be reached by the
328      !!         iterative solver
329      !!        c- apply the solver by a call to sol... routine
330      !!       -2- compute and add the free surface pressure gradient inclu-
331      !!      ding the additional force used to stabilize the equation.
332      !!
333      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
334      !!
335      !! References : Roullet and Madec 1999, JGR.
336      !!---------------------------------------------------------------------
337      !! * Arguments
338      INTEGER, INTENT( IN  ) :: &
339         &  kt         ! ocean time-step index
340      INTEGER, INTENT( OUT ) :: &
341         &  kindic     ! solver convergence flag (<0 if not converge)
342
343      !! * Local declarations
344      INTEGER  :: &
345         & ji, &     ! dummy loop indices
346         & jj, &
347         & jk
348      REAL(wp) :: &
349         & z2dt,   & ! temporary scalars
350         & z2dtg,  &
351         & zgcb,   &
352         & zbtd,   &
353         & ztdgu,  &
354         & ztdgv
355      !!----------------------------------------------------------------------
356      !
357      !
358      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt_adj')
359      !
360      IF( kt == nitend ) THEN
361         IF(lwp) WRITE(numout,*)
362         IF(lwp) WRITE(numout,*) 'dyn_spg_flt_adj : surface pressure gradient trend'
363         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   (free surface constant volume case)'
364      ENDIF
365
366      ! Local constant initialization
367      IF     ( neuler == 0 .AND. kt == nit000     ) THEN
368
369         z2dt = rdt             ! time step: Euler if restart from rest
370         CALL sol_mat(kt)       ! initialize matrix
371
372      ELSEIF (   neuler == 0 .AND. kt == nit000 + 1 ) THEN
373
374         z2dt = 2.0_wp * rdt    ! time step: leap-frog
375         CALL sol_mat(kt)       ! reinitialize matrix
376
377      ELSEIF (                   kt == nitend     ) THEN
378
379         z2dt = 2.0_wp * rdt    ! time step: leap-frog
380         CALL sol_mat(kt)       ! reinitialize matrix
381
382      ELSEIF ( neuler /= 0 .AND. kt == nit000     ) THEN
383
384         z2dt = 2.0_wp * rdt    ! time step: leap-frog
385         CALL sol_mat(kt)       ! initialize matrix
386
387      ELSE
388
389         z2dt = 2.0_wp * rdt    ! time step: leap-frog
390
391      ENDIF
392
393      z2dtg  = grav * z2dt
394
395      ! Add the trends multiplied by z2dt to the after velocity
396      ! -----------------------------------------------------------
397      !     ( c a u t i o n : (ua_ad,va_ad) here are the after velocity not the
398      !                       trend, the leap-frog time stepping will not
399      !                       be done in dynnxt_adj.F90 routine)
400      DO jk = 1, jpkm1
401         DO jj = 2, jpjm1
402            DO ji = fs_2, fs_jpim1   ! vector opt.
403               ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk)
404               va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk)
405               spgu_ad(ji,jj)  = spgu_ad(ji,jj)  + ua_ad(ji,jj,jk)
406               spgv_ad(ji,jj)  = spgv_ad(ji,jj)  + va_ad(ji,jj,jk)
407            END DO
408         END DO
409      END DO
410
411      ! Transport divergence gradient multiplied by z2dt
412      ! --------------------------------------------====
413      DO jj = 2, jpjm1
414         DO ji = fs_2, fs_jpim1   ! vector opt.
415            ! multiplied by z2dt
416            ztdgu = z2dt * spgu_ad(ji,jj)
417            ztdgv = z2dt * spgv_ad(ji,jj)
418            spgu_ad(ji,jj) = 0.0_wp
419            spgv_ad(ji,jj) = 0.0_wp
420            ! trend of Transport divergence gradient
421            ztdgu = ztdgu * z2dtg / e1u(ji,jj)
422            ztdgv = ztdgv * z2dtg / e2v(ji,jj)
423            gcx_ad(ji  ,jj  ) = gcx_ad(ji  ,jj  ) - ztdgu - ztdgv
424            gcx_ad(ji  ,jj+1) = gcx_ad(ji  ,jj+1) + ztdgv
425            gcx_ad(ji+1,jj  ) = gcx_ad(ji+1,jj  ) + ztdgu
426         END DO
427      END DO
428
429      ! Evaluate the next transport divergence
430      ! --------------------------------------
431      !    Iterative solver for the elliptic equation (except IF sol.=0)
432      !    (output in gcx_ad with boundary conditions applied)
433
434      kindic = 0
435      ncut = 0    !  Force solver
436      IF( ncut == 0 ) THEN
437         IF    ( nn_solv == 1 ) THEN ;  CALL ctl_stop('sol_pcg_adj not available in TMA yet')   ! diagonal preconditioned conjuguate gradient
438         ELSEIF( nn_solv == 2 ) THEN ;  CALL sol_sor_adj( kt, kindic )   ! successive-over-relaxation
439         ENDIF
440      ENDIF
441
442      ! Right hand side of the elliptic equation and first guess
443      ! --------------------------------------------------------
444      ! apply the lateral boundary conditions
445      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) &
446         &    CALL lbc_lnk_e_adj( gcb_ad, c_solver_pt, 1.0_wp )
447
448      DO jj = 2, jpjm1
449         DO ji = fs_2, fs_jpim1   ! vector opt.
450            ! First guess of the after barotropic transport divergence
451            zbtd = gcxb_ad(ji,jj) + 2.0_wp * gcx_ad(ji,jj)
452            gcxb_ad(ji,jj) = - gcx_ad(ji,jj)
453            gcx_ad (ji,jj) = zbtd
454            ! Divergence of the after vertically averaged velocity
455            zgcb = gcb_ad(ji,jj) * gcdprc(ji,jj)
456            gcb_ad(ji,jj) = 0.0_wp
457            spgu_ad(ji-1,jj  ) = spgu_ad(ji-1,jj  ) - zgcb
458            spgu_ad(ji  ,jj  ) = spgu_ad(ji  ,jj  ) + zgcb
459            spgv_ad(ji  ,jj-1) = spgv_ad(ji  ,jj-1) - zgcb
460            spgv_ad(ji  ,jj  ) = spgv_ad(ji  ,jj  ) + zgcb
461         END DO
462      END DO
463
464      IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_adj: lk_vvl is not available' )
465
466      ! Boundary conditions on (spgu_ad,spgv_ad)
467      CALL lbc_lnk_adj( spgu_ad, 'U', -1.0_wp )
468      CALL lbc_lnk_adj( spgv_ad, 'V', -1.0_wp )
469
470      ! transport: multiplied by the horizontal scale factor
471      DO jj = 2,jpjm1
472         DO ji = fs_2,fs_jpim1   ! vector opt.
473            spgu_ad(ji,jj) = spgu_ad(ji,jj) * e2u(ji,jj)
474            spgv_ad(ji,jj) = spgv_ad(ji,jj) * e1v(ji,jj)
475         END DO
476      END DO
477
478      ! compute the next vertically averaged velocity (effect of the additional force not included)
479      ! ---------------------------------------------
480
481      ! vertical sum
482!CDIR NOLOOPCHG
483      IF( lk_vopt_loop ) THEN     ! vector opt., forced unroll
484         DO jk = 1, jpkm1
485            DO ji = 1, jpij
486               ua_ad(ji,1,jk) = ua_ad(ji,1,jk) + fse3u(ji,1,jk) * spgu_ad(ji,1)
487               va_ad(ji,1,jk) = va_ad(ji,1,jk) + fse3v(ji,1,jk) * spgv_ad(ji,1)
488            END DO
489         END DO
490      ELSE                        ! No  vector opt.
491         DO jk = 1, jpkm1
492            DO jj = 2, jpjm1
493               DO ji = 2, jpim1
494                  ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + fse3u(ji,jj,jk) * spgu_ad(ji,jj)
495                  va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + fse3v(ji,jj,jk) * spgv_ad(ji,jj)
496               END DO
497            END DO
498         END DO
499      ENDIF
500
501      DO jj = 2, jpjm1
502         DO ji = fs_2, fs_jpim1 ! vector opt.
503            spgu_ad(ji,jj) = 0.0_wp
504            spgv_ad(ji,jj) = 0.0_wp
505         END DO
506      END DO
507
508      IF( nn_cla == 1 )   CALL cla_dynspg_adj( kt )      ! Cross Land Advection (update (ua,va))
509
510      ! Evaluate the masked next velocity (effect of the additional force not included)
511      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
512         !
513         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
514            DO jk = 1, jpkm1
515               DO jj = 2, jpjm1
516                  DO ji = fs_2, fs_jpim1   ! vector opt.
517                     ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + z2dt * ua_ad(ji,jj,jk) * umask(ji,jj,jk)
518                     ua_ad(ji,jj,jk) = z2dt * ua_ad(ji,jj,jk) * umask(ji,jj,jk)
519                     vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + z2dt * va_ad(ji,jj,jk) * vmask(ji,jj,jk)
520                     va_ad(ji,jj,jk) = z2dt * va_ad(ji,jj,jk) * vmask(ji,jj,jk)
521                  END DO
522               END DO
523            END DO
524            !
525         ELSE                          ! flux form : applied on thickness weighted velocity
526            DO jk = 1, jpkm1
527               DO jj = 2, jpjm1
528                  DO ji = fs_2, fs_jpim1   ! vector opt.
529                     ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
530                     ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk) * fse3u_b(ji,jj,jk)
531                     ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * z2dt * fse3u_n(ji,jj,jk)
532                     va_ad(ji,jj,jk) = va_ad(ji,jj,jk) / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
533                     vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk) * fse3v_b(ji,jj,jk)
534                     va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * z2dt * fse3v_n(ji,jj,jk)
535                 END DO
536               END DO
537            END DO
538            !
539         ENDIF
540         !
541      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
542         !
543         DO jk = 1, jpkm1               ! unweighted time stepping
544            DO jj = 2, jpjm1
545               DO ji = fs_2, fs_jpim1   ! vector opt.
546                  ua_ad(  ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk)
547                  ub_ad(  ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk)
548                  spgu_ad(ji,jj   ) = spgu_ad(ji,jj)  + ua_ad(ji,jj,jk) * z2dt
549                  ua_ad(  ji,jj,jk) = ua_ad(ji,jj,jk) * z2dt
550                  va_ad(  ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk)
551                  vb_ad(  ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk)
552                  spgv_ad(ji,jj   ) = spgv_ad(ji,jj)  + va_ad(ji,jj,jk) * z2dt
553                  va_ad(  ji,jj,jk) = va_ad(ji,jj,jk) * z2dt
554               END DO
555            END DO
556         END DO
557         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
558            DO ji = fs_2, fs_jpim1     ! vector opt.
559               spgu_ad(ji  ,jj  ) = spgu_ad(ji  ,jj  ) * grav / e1u(ji,jj)
560               spgv_ad(ji  ,jj  ) = spgv_ad(ji  ,jj  ) * grav / e2v(ji,jj)
561               sshn_ad(ji  ,jj  ) = sshn_ad(ji  ,jj  ) + spgv_ad(ji,jj)
562               sshn_ad(ji  ,jj+1) = sshn_ad(ji  ,jj+1) - spgv_ad(ji,jj)
563               sshn_ad(ji  ,jj  ) = sshn_ad(ji  ,jj  ) + spgu_ad(ji,jj)
564               sshn_ad(ji+1,jj  ) = sshn_ad(ji+1,jj  ) - spgu_ad(ji,jj)
565               spgu_ad(ji  ,jj  ) = 0.0_wp
566               spgv_ad(ji  ,jj  ) = 0.0_wp
567            END DO
568         END DO
569      ENDIF
570
571      IF( kt == nit000 ) THEN
572         ! set to zero free surface specific arrays
573         spgu_ad(:,:) = 0.0_wp                    ! surface pressure gradient (i-direction)
574         spgv_ad(:,:) = 0.0_wp                    ! surface pressure gradient (j-direction)
575      ENDIF
576      !
577      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt_adj')
578      !
579   END SUBROUTINE dyn_spg_flt_adj
580
581   SUBROUTINE dyn_spg_flt_adj_tst( kumadt )
582      !!-----------------------------------------------------------------------
583      !!
584      !!                  ***  ROUTINE dyn_spg_flt_adj_tst ***
585      !!
586      !! ** Purpose : Test the adjoint routine.
587      !!
588      !! ** Method  : Verify the scalar product
589      !!
590      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
591      !!
592      !!              where  L   = tangent routine
593      !!                     L^T = adjoint routine
594      !!                     W   = diagonal matrix of scale factors
595      !!                     dx  = input perturbation (random field)
596      !!                     dy  = L dx
597      !!
598      !! ** Action  :
599      !!
600      !! History :
601      !!        ! 09-01 (A. Weaver)
602      !!-----------------------------------------------------------------------
603      !! * Modules used
604
605      !! * Arguments
606      INTEGER, INTENT(IN) :: &
607         & kumadt        ! Output unit
608
609      !! * Local declarations
610      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
611         & zua_tlin,    & ! Tangent input: ua_tl
612         & zva_tlin,    & ! Tangent input: va_tl
613         & zub_tlin,    & ! Tangent input: ub_tl
614         & zvb_tlin,    & ! Tangent input: vb_tl
615         & zua_tlout,   & ! Tangent output: ua_tl
616         & zva_tlout,   & ! Tangent output: va_tl
617         & zub_tlout,   & ! Tangent output: ua_tl
618         & zvb_tlout,   & ! Tangent output: va_tl
619         & zub_adin,    & ! Tangent output: ua_ad
620         & zvb_adin,    & ! Tangent output: va_ad
621         & zua_adin,    & ! Adjoint input: ua_ad
622         & zva_adin,    & ! Adjoint input: va_ad
623         & zua_adout,   & ! Adjoint output: ua_ad
624         & zva_adout,   & ! Adjoint output: va_ad
625         & zub_adout,   & ! Adjoint oputput: ub_ad
626         & zvb_adout,   & ! Adjoint output: vb_ad
627         & znu            ! 3D random field for u
628
629      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
630         & zgcx_tlin, zgcxb_tlin, zgcb_tlin, zgcx_tlout, zgcxb_tlout, zgcb_tlout,  &
631         & zgcx_adin, zgcxb_adin, zgcb_adin, zgcx_adout, zgcxb_adout, zgcb_adout,  &
632         & zspgu_tlout, zspgv_tlout, zspgu_adin, zspgv_adin
633
634      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
635         & zsshn_tlin, & ! Tangent input: sshn_tl
636         & zsshn_adout,& ! Adjoint output: sshn_ad
637         & zemp_tlin,  & ! Tangent input: emp_tl
638         & zemp_adout, & ! Adjoint output: emp_ad
639         & znssh         ! 2D random field for SSH
640      REAL(wp) :: &
641         & zsp1,    &   ! scalar product involving the tangent routine
642         & zsp2         ! scalar product involving the adjoint routine
643      INTEGER :: &
644         & indic, &
645         & istp
646      INTEGER :: &
647         & ji, &
648         & jj, &
649         & jk, &
650         & kmod, &
651         & jstp
652      CHARACTER (LEN=14) :: &
653         & cl_name
654      INTEGER ::             &
655         & jpert
656      INTEGER, PARAMETER :: jpertmax = 7
657
658      ! Allocate memory
659
660      ALLOCATE( &
661         & zua_tlin(jpi,jpj,jpk),  &
662         & zva_tlin(jpi,jpj,jpk),  &
663         & zub_tlin(jpi,jpj,jpk),  &
664         & zvb_tlin(jpi,jpj,jpk),  &
665         & zua_tlout(jpi,jpj,jpk), &
666         & zva_tlout(jpi,jpj,jpk), &
667         & zua_adin(jpi,jpj,jpk),  &
668         & zva_adin(jpi,jpj,jpk),  &
669         & zua_adout(jpi,jpj,jpk), &
670         & zva_adout(jpi,jpj,jpk), &
671         & zub_adout(jpi,jpj,jpk), &
672         & zvb_adout(jpi,jpj,jpk), &
673         & znu(jpi,jpj,jpk)        &
674         & )
675      ALLOCATE( &
676         & zsshn_tlin(jpi,jpj), &
677         & zsshn_adout(jpi,jpj),&
678         & zemp_tlin(jpi,jpj),  &
679         & zemp_adout(jpi,jpj), &
680         & znssh(jpi,jpj)       &
681         & )
682
683      ALLOCATE( zgcx_tlin (jpi,jpj), zgcx_tlout (jpi,jpj), zgcx_adin (jpi,jpj), zgcx_adout (jpi,jpj),  &
684                zgcxb_tlin(jpi,jpj), zgcxb_tlout(jpi,jpj), zgcxb_adin(jpi,jpj), zgcxb_adout(jpi,jpj),  &
685                zgcb_tlin (jpi,jpj), zgcb_tlout (jpi,jpj), zgcb_adin (jpi,jpj), zgcb_adout (jpi,jpj)   &
686               & )
687
688      ALLOCATE ( zub_tlout(jpi,jpj,jpk), zvb_tlout(jpi,jpj,jpk),  &
689                 zub_adin (jpi,jpj,jpk), zvb_adin (jpi,jpj,jpk)  )
690
691      ALLOCATE( zspgu_tlout (jpi,jpj), zspgv_tlout (jpi,jpj), zspgu_adin (jpi,jpj), zspgv_adin (jpi,jpj))
692
693      !=========================================================================
694      !     dx = ( sshb_tl, sshn_tl, ub_tl, ua_tl, vb_tl, va_tl, wn_tl, emp_tl )
695      ! and dy = ( sshb_tl, sshn_tl, ua_tl, va_tl )
696      !=========================================================================
697
698      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
699
700      DO jstp = nit000, nit000 + 1
701         DO jpert = 1, jpertmax
702            istp = jstp
703
704            !--------------------------------------------------------------------
705            ! Reset the tangent and adjoint variables
706            !--------------------------------------------------------------------
707
708            zub_tlin (:,:,:) = 0.0_wp
709            zvb_tlin (:,:,:) = 0.0_wp
710            zua_tlin (:,:,:) = 0.0_wp
711            zva_tlin (:,:,:) = 0.0_wp
712            zua_tlout(:,:,:) = 0.0_wp
713            zva_tlout(:,:,:) = 0.0_wp
714            zua_adin (:,:,:) = 0.0_wp
715            zva_adin (:,:,:) = 0.0_wp
716            zub_adout(:,:,:) = 0.0_wp
717            zvb_adout(:,:,:) = 0.0_wp
718            zua_adout(:,:,:) = 0.0_wp
719            zva_adout(:,:,:) = 0.0_wp
720
721            zsshn_tlin (:,:) = 0.0_wp
722            zemp_tlin  (:,:) = 0.0_wp
723            zsshn_adout(:,:) = 0.0_wp
724            zemp_adout (:,:) = 0.0_wp
725            zspgu_adin (:,:) = 0.0_wp
726            zspgv_adin (:,:) = 0.0_wp
727            zspgu_tlout(:,:) = 0.0_wp
728            zspgv_tlout(:,:) = 0.0_wp
729
730            zgcx_tlout (:,:) = 0.0_wp ; zgcx_adin (:,:) = 0.0_wp ; zgcx_adout (:,:) = 0.0_wp
731            zgcxb_tlout(:,:) = 0.0_wp ; zgcxb_adin(:,:) = 0.0_wp ; zgcxb_adout(:,:) = 0.0_wp
732            zgcb_tlout (:,:) = 0.0_wp ; zgcb_adin (:,:) = 0.0_wp ; zgcb_adout (:,:) = 0.0_wp
733
734            ub_tl(:,:,:) = 0.0_wp
735            vb_tl(:,:,:) = 0.0_wp
736            ua_tl(:,:,:) = 0.0_wp
737            va_tl(:,:,:) = 0.0_wp
738            sshn_tl(:,:) = 0.0_wp
739            emp_tl(:,:)  = 0.0_wp
740            gcb_tl(:,:)  = 0.0_wp
741            gcx_tl(:,:)  = 0.0_wp
742            gcxb_tl(:,:) = 0.0_wp
743            spgu_tl(:,:) = 0.0_wp
744            spgv_tl(:,:) = 0.0_wp
745            ub_ad(:,:,:) = 0.0_wp
746            vb_ad(:,:,:) = 0.0_wp
747            ua_ad(:,:,:) = 0.0_wp
748            va_ad(:,:,:) = 0.0_wp
749            sshn_ad(:,:) = 0.0_wp
750            emp_ad(:,:)  = 0.0_wp
751            gcb_ad(:,:)  = 0.0_wp
752            gcx_ad(:,:)  = 0.0_wp
753            gcxb_ad(:,:) = 0.0_wp
754            spgu_ad(:,:) = 0.0_wp
755            spgv_ad(:,:) = 0.0_wp
756            !--------------------------------------------------------------------
757            ! Initialize the tangent input with random noise: dx
758            !--------------------------------------------------------------------
759            IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN
760
761               CALL grid_random(  znu, 'U', 0.0_wp, stdu )
762
763               DO jk = 1, jpk
764                  DO jj = nldj, nlej
765                     DO ji = nldi, nlei
766                        zua_tlin(ji,jj,jk) = znu(ji,jj,jk)
767                     END DO
768                  END DO
769               END DO
770
771            ENDIF
772            IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN
773               CALL grid_random(  znu, 'V', 0.0_wp, stdv )
774
775               DO jk = 1, jpk
776                  DO jj = nldj, nlej
777                     DO ji = nldi, nlei
778                        zva_tlin(ji,jj,jk) = znu(ji,jj,jk)
779                     END DO
780                  END DO
781               END DO
782
783            ENDIF
784            IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN
785               CALL grid_random(  znu, 'U', 0.0_wp, stdu )
786
787               DO jk = 1, jpk
788                  DO jj = nldj, nlej
789                     DO ji = nldi, nlei
790                        zub_tlin(ji,jj,jk) = znu(ji,jj,jk)
791                     END DO
792                  END DO
793               END DO
794
795            ENDIF
796            IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN
797               CALL grid_random(  znu, 'V', 0.0_wp, stdv )
798
799               DO jk = 1, jpk
800                  DO jj = nldj, nlej
801                     DO ji = nldi, nlei
802                        zvb_tlin(ji,jj,jk) = znu(ji,jj,jk)
803                     END DO
804                  END DO
805               END DO
806
807            ENDIF
808            IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN
809               CALL grid_random(  znssh, 'T', 0.0_wp, stdemp )
810
811               DO jj = nldj, nlej
812                  DO ji = nldi, nlei
813                     zemp_tlin(ji,jj) = znssh(ji,jj)
814                  END DO
815               END DO
816
817            ENDIF
818            IF ( (jpert == 6) .OR. (jpert == jpertmax) ) THEN
819               CALL grid_random(  znssh, 'T', 0.0_wp, stdssh )
820               DO jj = nldj, nlej
821                  DO ji = nldi, nlei
822                     zsshn_tlin(ji,jj) = znssh(ji,jj)
823                  END DO
824               END DO
825            END IF
826            zgcx_tlin  (:,:) = ( zua_tlin(:,:,1) + zub_tlin(:,:,1) ) / 10.
827            zgcxb_tlin (:,:) = ( zua_tlin(:,:,2) + zub_tlin(:,:,2) ) / 10.
828            !--------------------------------------------------------------------
829            ! Call the tangent routine: dy = L dx
830            !--------------------------------------------------------------------
831
832            ua_tl(:,:,:) = zua_tlin(:,:,:)
833            va_tl(:,:,:) = zva_tlin(:,:,:)
834            ub_tl(:,:,:) = zub_tlin(:,:,:)
835            vb_tl(:,:,:) = zvb_tlin(:,:,:)
836            emp_tl (:,:) = zemp_tlin (:,:)
837            sshn_tl(:,:) = zsshn_tlin(:,:)
838
839            gcx_tl (:,:) = 0.e0              ;   gcxb_tl(:,:) = 0.e0
840            gcb_tl (:,:) = 0.e0
841            gcx_tl (:,:) = zgcx_tlin (:,:)   ;   gcxb_tl(:,:) = zgcxb_tlin(:,:)
842
843            CALL dyn_spg_flt_tan( istp, indic )
844
845            zua_tlout(:,:,:) = ua_tl(:,:,:)   ;   zva_tlout(:,:,:) = va_tl(:,:,:)
846            zspgu_tlout(:,:) = spgu_tl(:,:)   ;   zspgv_tlout(:,:) = spgv_tl(:,:)
847            zgcb_tlout (:,:) = gcb_tl (:,:)
848
849            !--------------------------------------------------------------------
850            ! Initialize the adjoint variables: dy^* = W dy
851            !--------------------------------------------------------------------
852
853            DO jk = 1, jpk
854               DO jj = nldj, nlej
855                  DO ji = nldi, nlei
856                     zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
857                        &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
858                        &               * umask(ji,jj,jk)
859                     zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
860                        &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
861                        &               * vmask(ji,jj,jk)
862                  END DO
863               END DO
864            END DO
865            DO jj = nldj, nlej
866               DO ji = nldi, nlei
867                  zgcb_adin (ji,jj) = zgcb_tlout (ji,jj) &
868                     &              * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,1) * tmask(ji,jj,1)
869                  zspgu_adin (ji,jj) = zspgu_tlout (ji,jj) &
870                     &              * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) * umask(ji,jj,1)
871                  zspgv_adin(ji,jj) = zspgv_tlout(ji,jj) &
872                     &              * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) * vmask(ji,jj,1)
873               END DO
874            END DO
875
876            !--------------------------------------------------------------------
877            ! Compute the scalar product: ( L dx )^T W dy
878            !--------------------------------------------------------------------
879
880            zsp1 =   DOT_PRODUCT( zua_tlout  , zua_adin   ) &
881               &   + DOT_PRODUCT( zgcb_tlout , zgcb_adin  ) &
882               &   + DOT_PRODUCT( zspgu_tlout , zspgu_adin  ) &
883               &   + DOT_PRODUCT( zspgv_tlout , zspgv_adin  ) &
884               &   + DOT_PRODUCT( zva_tlout  , zva_adin   )
885
886
887            !--------------------------------------------------------------------
888            ! Call the adjoint routine: dx^* = L^T dy^*
889            !--------------------------------------------------------------------
890
891            ua_ad(:,:,:) = zua_adin(:,:,:)
892            va_ad(:,:,:) = zva_adin(:,:,:)
893
894            gcx_ad (:,:)   = 0.0_wp            ;   gcxb_ad(:,:) = 0.0_wp
895            gcb_ad (:,:)   = zgcb_adin (:,:)
896            spgu_ad(:,:)   = zspgu_adin(:,:)
897            spgv_ad(:,:)   = zspgv_adin(:,:)
898            ub_ad  (:,:,:) = zub_adin  (:,:,:) ;  vb_ad  (:,:,:) = zvb_adin  (:,:,:)
899
900            CALL dyn_spg_flt_adj( istp, indic )
901
902            zua_adout(:,:,:) = ua_ad(:,:,:)
903            zva_adout(:,:,:) = va_ad(:,:,:)
904            zub_adout(:,:,:) = ub_ad(:,:,:)
905            zvb_adout(:,:,:) = vb_ad(:,:,:)
906            zemp_adout (:,:) = emp_ad (:,:)
907            zsshn_adout(:,:) = sshn_ad(:,:)
908            zgcx_adout (:,:) = gcx_ad (:,:)
909            zgcxb_adout(:,:) = gcxb_ad(:,:)
910
911            !--------------------------------------------------------------------
912            ! Compute the scalar product: dx^T L^T W dy
913            !--------------------------------------------------------------------
914
915            zsp2 =   DOT_PRODUCT( zua_tlin  , zua_adout   ) &
916               &   + DOT_PRODUCT( zva_tlin  , zva_adout   ) &
917               &   + DOT_PRODUCT( zub_tlin  , zub_adout   ) &
918               &   + DOT_PRODUCT( zvb_tlin  , zvb_adout   ) &
919               &   + DOT_PRODUCT( zgcx_tlin , zgcx_adout  ) &
920               &   + DOT_PRODUCT( zgcxb_tlin, zgcxb_adout ) &
921               &   + DOT_PRODUCT( zsshn_tlin, zsshn_adout ) &
922               &   + DOT_PRODUCT( zemp_tlin , zemp_adout  )
923
924            ! Compare the scalar products
925
926            !    14 char:'12345678901234'
927            IF ( istp == nit000 ) THEN
928               SELECT CASE (jpert)
929               CASE(1)
930                  cl_name = 'spg_flt  Ua T1'
931               CASE(2)
932                  cl_name = 'spg_flt  Va T1'
933               CASE(3)
934                  cl_name = 'spg_flt  Ub T1'
935               CASE(4)
936                  cl_name = 'spg_flt  Vb T1'
937               CASE(5)
938                  cl_name = 'spg_flt emp T1'
939               CASE(6)
940                  cl_name = 'spg_flt ssh T1'
941               CASE(jpertmax)
942                  cl_name = 'dyn_spg_flt T1'
943               END SELECT
944            ELSEIF ( istp == nit000 + 1 ) THEN
945               SELECT CASE (jpert)
946               CASE(1)
947                  cl_name = 'spg_flt  Ua T2'
948               CASE(2)
949                  cl_name = 'spg_flt  Va T2'
950               CASE(3)
951                  cl_name = 'spg_flt  Ub T2'
952               CASE(4)
953                  cl_name = 'spg_flt  Vb T2'
954               CASE(5)
955                  cl_name = 'spg_flt emp T2'
956               CASE(6)
957                  cl_name = 'spg_flt ssh T2'
958               CASE(jpertmax)
959                  cl_name = 'dyn_spg_flt T2'
960               END SELECT
961            END IF
962            CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
963
964         END DO
965      END DO
966
967!nn_nmod = kmod  ! restore initial frequency of test for the SOR solver
968
969      ! Deallocate memory
970
971      DEALLOCATE( &
972         & zua_tlin,  &
973         & zva_tlin,  &
974         & zub_tlin,  &
975         & zvb_tlin,  &
976         & zua_tlout, &
977         & zva_tlout, &
978         & zua_adin,  &
979         & zva_adin,  &
980         & zua_adout, &
981         & zva_adout, &
982         & zub_adout, &
983         & zvb_adout, &
984         & znu        &
985         & )
986      DEALLOCATE( &
987         & zsshn_tlin, &
988         & zemp_tlin,  &
989         & zsshn_adout,&
990         & zemp_adout, &
991         & znssh       &
992         & )
993      DEALLOCATE( zgcx_tlin , zgcx_tlout , zgcx_adin , zgcx_adout,  &
994         & zgcxb_tlin, zgcxb_tlout, zgcxb_adin, zgcxb_adout,  &
995         & zgcb_tlin , zgcb_tlout , zgcb_adin , zgcb_adout    &
996         & )
997      DEALLOCATE ( zub_tlout, zvb_tlout, zub_adin , zvb_adin )
998   END SUBROUTINE dyn_spg_flt_adj_tst
999
1000# else
1001   !!----------------------------------------------------------------------
1002   !!   Default case :   Empty module   No standart explicit free surface
1003   !!----------------------------------------------------------------------
1004CONTAINS
1005   SUBROUTINE dyn_spg_flt_tan( kt, kindic )       ! Empty routine
1006      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt
1007   END SUBROUTINE dyn_spg_flt_tan
1008   SUBROUTINE dyn_spg_flt_adj( kt, kindic )       ! Empty routine
1009      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt
1010   END SUBROUTINE dyn_spg_flt_adj
1011   SUBROUTINE dyn_spg_flt_adj_tst( kt )       ! Empty routine
1012      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt
1013   END SUBROUTINE dyn_spg_flt_adj_tst
1014   SUBROUTINE dyn_spg_flt_tlm_tst( kt )       ! Empty routine
1015      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt
1016   END SUBROUTINE dyn_spg_flt_tlm_tst
1017# endif
1018#endif
1019END MODULE dynspg_flt_tam
Note: See TracBrowser for help on using the repository browser.