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_exp_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_exp_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: 13.8 KB
Line 
1MODULE dynspg_exp_tam
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_exp_tam  TANGENT/ADJOINT OF MODULE dynspg_exp***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6   !! History of the direct module:
7   !!            2.0  !  2005-11  (V. Garnier, G. Madec, L. Bessieres) Original code
8   !!            3.2  !  2009-06  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
9   !! History of the tam module:
10   !!            3.2  !  2010-06  (A. Vidard) tam of the 2009-06 version
11   !!----------------------------------------------------------------------
12#if defined key_dynspg_exp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_dynspg_exp'                              explicit free surface
15   !!----------------------------------------------------------------------
16   !!   dyn_spg_exp  : update the momentum trend with the surface
17   !!                      pressure gradient in the free surface constant
18   !!                      volume case with vector optimization
19   !!----------------------------------------------------------------------
20   USE par_kind
21   USE phycst
22   USE par_oce
23   USE oce_tam
24   USE dom_oce
25   USE gridrandom
26   USE dotprodfld
27   USE paresp
28   USE in_out_manager
29   USE tstool_tam
30   USE timing          ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_spg_exp_tan       ! routine called by step.F90
36   PUBLIC   dyn_spg_exp_adj       ! routine called by step.F90
37   PUBLIC   dyn_spg_exp_adj_tst   ! routine called by tamtst.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
44   !! $Id$
45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE dyn_spg_exp_tan( kt )
51      !!----------------------------------------------------------------------
52      !!                  ***  routine dyn_spg_exp_tan  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the surface pressure
55      !!              gradient in case of explicit free surface formulation and
56      !!              add it to the general trend of momentum equation.
57      !!
58      !! ** Method  :   Explicit free surface formulation. Add to the general
59      !!              momentum trend the surface pressure gradient :
60      !!                      (ua,va) = (ua,va) + (spgu,spgv)
61      !!              where spgu = -1/rau0 d/dx(ps) = -g/e1u di( sshn )
62      !!                    spgv = -1/rau0 d/dy(ps) = -g/e2v dj( sshn )
63      !!
64      !! ** Action :   (ua,va)   trend of horizontal velocity increased by
65      !!                         the surf. pressure gradient trend
66      !!---------------------------------------------------------------------
67      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
68      !!
69      INTEGER  ::   ji, jj, jk               ! dummy loop indices
70      !!----------------------------------------------------------------------
71      !
72      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_exp_tan')
73      !
74      IF( kt == nit000 ) THEN
75         IF(lwp) WRITE(numout,*)
76         IF(lwp) WRITE(numout,*) 'dyn_spg_exp_tan : surface pressure gradient trend'
77         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (explicit free surface)'
78         !
79         spgu_tl(:,:) = 0._wp  ;   spgv_tl(:,:) = 0._wp
80         !
81         IF( lk_vvl .AND. lwp ) WRITE(numout,*) '              lk_vvl=T : spg is included in dynhpg'
82      ENDIF
83      IF( .NOT. lk_vvl ) THEN          !* fixed volume : add the surface pressure gradient trend
84         !
85         DO jj = 2, jpjm1                    ! now surface pressure gradient
86            DO ji = fs_2, fs_jpim1   ! vector opt.
87               spgu_tl(ji,jj) = - grav * ( sshn_tl(ji+1,jj) - sshn_tl(ji,jj) ) / e1u(ji,jj)
88               spgv_tl(ji,jj) = - grav * ( sshn_tl(ji,jj+1) - sshn_tl(ji,jj) ) / e2v(ji,jj)
89            END DO
90         END DO
91         DO jk = 1, jpkm1                    ! Add it to the general trend
92            DO jj = 2, jpjm1
93               DO ji = fs_2, fs_jpim1   ! vector opt.
94                  ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + spgu_tl(ji,jj)
95                  va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + spgv_tl(ji,jj)
96               END DO
97            END DO
98         END DO
99         !
100      ENDIF
101      !
102      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_exp_tan')
103      !
104   END SUBROUTINE dyn_spg_exp_tan
105   SUBROUTINE dyn_spg_exp_adj( kt )
106      !!----------------------------------------------------------------------
107      !!                  ***  routine dyn_spg_exp_adj  ***
108      !!
109      !! ** Purpose :   Compute the now trend due to the surface pressure
110      !!              gradient in case of explicit free surface formulation and
111      !!              add it to the general trend of momentum equation.
112      !!
113      !! ** Method  :   Explicit free surface formulation. Add to the general
114      !!              momentum trend the surface pressure gradient :
115      !!                      (ua,va) = (ua,va) + (spgu,spgv)
116      !!              where spgu = -1/rau0 d/dx(ps) = -g/e1u di( sshn )
117      !!                    spgv = -1/rau0 d/dy(ps) = -g/e2v dj( sshn )
118      !!
119      !! ** Action :   (ua,va)   trend of horizontal velocity increased by
120      !!                         the surf. pressure gradient trend
121      !!---------------------------------------------------------------------
122      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
123      !!
124      INTEGER  ::   ji, jj, jk               ! dummy loop indices
125      !!----------------------------------------------------------------------
126      !
127      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_exp_adj')
128      !
129      IF( kt == nitend ) THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) 'dyn_spg_exp_adj : surface pressure gradient trend'
132         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (explicit free surface)'
133      END IF
134
135      spgu_ad(:,:) = 0._wp   ;   spgv_ad(:,:) = 0._wp
136
137      IF( .NOT. lk_vvl ) THEN          !* fixed volume : add the surface pressure gradient trend
138         !
139         DO jk = 1, jpkm1                    ! Add it to the general trend
140            DO jj = 2, jpjm1
141               DO ji = fs_2, fs_jpim1   ! vector opt.
142                  spgu_ad(ji,jj) = spgu_ad(ji,jj) + ua_ad(ji,jj,jk)
143                  spgv_ad(ji,jj) = spgv_ad(ji,jj) + va_ad(ji,jj,jk)
144               END DO
145            END DO
146         END DO
147         DO jj = jpjm1, 2, -1                    ! now surface pressure gradient
148            DO ji = fs_jpim1, fs_2, -1           ! vector opt.
149
150               spgu_ad(ji,jj) = - grav * spgu_ad(ji,jj) / e1u(ji,jj)
151               spgv_ad(ji,jj) = - grav * spgv_ad(ji,jj) / e2v(ji,jj)
152               sshn_ad(ji+1,jj) = sshn_ad(ji+1,jj) + spgu_ad(ji,jj)
153               sshn_ad(ji,jj+1) = sshn_ad(ji,jj+1) + spgv_ad(ji,jj)
154               sshn_ad(ji,jj)   = sshn_ad(ji,jj) - spgu_ad(ji,jj) - spgv_ad(ji,jj)
155            END DO
156         END DO
157         !
158      ENDIF
159      !
160   !
161   IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_exp_adj')
162   !
163   END SUBROUTINE dyn_spg_exp_adj
164   SUBROUTINE dyn_spg_exp_adj_tst( kumadt )
165      !!-----------------------------------------------------------------------
166      !!
167      !!                  ***  ROUTINE dyn_spg_exp_adj_tst ***
168      !!
169      !! ** Purpose : Test the adjoint routine.
170      !!
171      !! ** Method  : Verify the scalar product
172      !!
173      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
174      !!
175      !!              where  L   = tangent routine
176      !!                     L^T = adjoint routine
177      !!                     W   = diagonal matrix of scale factors
178      !!                     dx  = input perturbation (random field)
179      !!                     dy  = L dx
180      !!
181      !!
182      !! History :
183      !!        ! 2010-06 (A. Vidard)
184      !!-----------------------------------------------------------------------
185      !! * Modules used
186
187      !! * Arguments
188      INTEGER, INTENT(IN) :: &
189         & kumadt             ! Output unit
190
191      !! * Local declarations
192      INTEGER ::  &
193         & ji,    &        ! dummy loop indices
194         & jj,    &
195         & jk
196      INTEGER, DIMENSION(jpi,jpj) :: &
197         & iseed_2d        ! 2D seed for the random number generator
198      REAL(KIND=wp) :: &
199         & zsp1,         & ! scalar product involving the tangent routine
200         & zsp2            ! scalar product involving the adjoint routine
201      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
202         & zua_tlin ,     & ! Tangent input
203         & zva_tlin ,     & ! Tangent input
204         & zua_tlout,     & ! Tangent output
205         & zva_tlout,     & ! Tangent output
206         & zua_adin ,     & ! Adjoint input
207         & zva_adin ,     & ! Adjoint input
208         & zua_adout,     & ! Adjoint output
209         & zva_adout,     & ! Adjoint output
210         & zr3d           ! 3D random field
211     REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
212        & zsshn_tlin,    &
213        & zsshn_adout,   &
214        & zr2d
215       CHARACTER(LEN=14) :: &
216         & cl_name
217      ! Allocate memory
218
219      ALLOCATE( &
220         & zua_tlin(   jpi,jpj,jpk),     &
221         & zva_tlin(   jpi,jpj,jpk),     &
222         & zsshn_tlin( jpi,jpj    ),     &
223         & zua_tlout(  jpi,jpj,jpk),     &
224         & zva_tlout(  jpi,jpj,jpk),     &
225         & zua_adin(   jpi,jpj,jpk),     &
226         & zva_adin(   jpi,jpj,jpk),     &
227         & zua_adout(  jpi,jpj,jpk),     &
228         & zva_adout(  jpi,jpj,jpk),     &
229         & zsshn_adout(jpi,jpj    ),     &
230         & zr3d(       jpi,jpj,jpk),     &
231         & zr2d(       jpi,jpj    )      &
232         & )
233      !==================================================================
234      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
235      !    dy = ( hdivb_tl, hdivn_tl )
236      !==================================================================
237
238      !--------------------------------------------------------------------
239      ! Reset the tangent and adjoint variables
240      !--------------------------------------------------------------------
241
242      ua_ad( :,:,:) = 0.0_wp
243      va_ad( :,:,:) = 0.0_wp
244      sshn_ad( :,:) = 0.0_wp
245      !--------------------------------------------------------------------
246      ! Initialize the tangent input with random noise: dx
247      !--------------------------------------------------------------------
248
249      CALL grid_random( zr3d, 'U', 0.0_wp, stdu )
250      zua_tlin(:,:,:) = zr3d(:,:,:)
251      CALL grid_random( zr3d, 'V', 0.0_wp, stdv )
252      zva_tlin(:,:,:) = zr3d(:,:,:)
253      CALL grid_random( zr2d, 'T', 0.0_wp, stdssh )
254      zsshn_tlin(:,:) = zr2d(:,:)
255
256
257      ua_tl = zua_tlin
258      va_tl = zva_tlin
259      sshn_tl = zsshn_tlin
260      CALL dyn_spg_exp_tan( nit000 )
261      zua_tlout = ua_tl
262      zva_tlout = va_tl
263      !--------------------------------------------------------------------
264      ! Initialize the adjoint variables: dy^* = W dy
265      !--------------------------------------------------------------------
266
267      DO jk = 1, jpk
268        DO jj = nldj, nlej
269           DO ji = nldi, nlei
270              zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
271                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
272                 &               * umask(ji,jj,jk)
273              zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
274                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
275                 &               * vmask(ji,jj,jk)
276            END DO
277         END DO
278      END DO
279      !--------------------------------------------------------------------
280      ! Compute the scalar product: ( L dx )^T W dy
281      !--------------------------------------------------------------------
282
283      zsp1 = DOT_PRODUCT( zua_tlout, zua_adin ) &
284         & + DOT_PRODUCT( zva_tlout, zva_adin )
285
286      !--------------------------------------------------------------------
287      ! Call the adjoint routine: dx^* = L^T dy^*
288      !--------------------------------------------------------------------
289
290      ua_ad = zua_adin
291      va_ad = zva_adin
292
293      CALL dyn_spg_exp_adj( nit000 )
294
295      zua_adout   = ua_ad
296      zva_adout   = va_ad
297      zsshn_adout = sshn_ad
298
299      zsp2 = DOT_PRODUCT( zua_tlin  , zua_adout   ) &
300         & + DOT_PRODUCT( zva_tlin  , zva_adout   ) &
301         & + DOT_PRODUCT( zsshn_tlin, zsshn_adout )
302
303      ! 14 char:'12345678901234'
304      cl_name = 'dyn_spg_exp   '
305      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
306
307      DEALLOCATE(         &
308         & zua_tlin,      &
309         & zva_tlin,      &
310         & zsshn_tlin,    &
311         & zua_tlout,     &
312         & zva_tlout,     &
313         & zua_adin,      &
314         & zva_adin,      &
315         & zua_adout,     &
316         & zva_adout,     &
317         & zsshn_adout,   &
318         & zr3d,          &
319         & zr2d           &
320         & )
321
322   END SUBROUTINE dyn_spg_exp_adj_tst
323
324#else
325   !!----------------------------------------------------------------------
326   !!   Default case :   Empty module   No standart explicit free surface
327   !!----------------------------------------------------------------------
328CONTAINS
329   SUBROUTINE dyn_spg_exp_tan( kt )       ! Empty routine
330      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
331   END SUBROUTINE dyn_spg_exp_tan
332   SUBROUTINE dyn_spg_exp_adj( kt )       ! Empty routine
333      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
334   END SUBROUTINE dyn_spg_exp_adj
335   SUBROUTINE dyn_spg_exp_adj_tst( kt )       ! Empty routine
336      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
337   END SUBROUTINE dyn_spg_exp_adj_tst
338#endif
339
340   !!======================================================================
341END MODULE dynspg_exp_tam
Note: See TracBrowser for help on using the repository browser.