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.
dynadv_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/dynadv_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: 17.7 KB
Line 
1MODULE dynadv_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                       ***  MODULE  dynadv_tam  ***
5   !! Ocean active tracers:  advection scheme control
6   !!                        Tangent and Adjoint module
7   !!==============================================================================
8   !! History of the direct module:
9   !!           9.0  ! 2006-11  (G. Madec)  Original code
10   !! History of the TAM module:
11   !!           9.0  ! 2008-08  (A. Vidard) first version
12   !!   NEMO    3.2  ! 2010-04 (F. Vigilant) 3.2 version
13   !!   NEMO    3.4  ! 2012-07 (P.-A. Bouttier) 3.4 version
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   dyn_adv      : compute the momentum advection trend
18   !!   dyn_adv_ctl  : control the different options of advection scheme
19   !!----------------------------------------------------------------------
20   USE par_kind
21   USE par_oce
22   USE oce
23   USE dom_oce
24   USE oce_tam
25   USE in_out_manager
26   USE gridrandom
27   USE dotprodfld
28   USE tstool_tam
29   USE dynadv
30   USE dynadv_cen2_tam
31   USE dynadv_ubs_tam
32   USE dynkeg_tam
33   USE dynzad_tam
34   USE sshwzv_tam
35   USE sshwzv
36   USE divcur
37   USE divcur_tam
38   USE in_out_manager
39   USE lib_mpp
40   USE timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC dyn_adv_tan     ! routine called by steptan module
46   PUBLIC dyn_adv_adj     ! routine called by stepadj module
47   PUBLIC dyn_adv_adj_tst ! routine called by the tst module
48   PUBLIC dyn_adv_init_tam
49#if defined key_tst_tlm
50   PUBLIC dyn_adv_tlm_tst ! routine called by tamtst.F90
51#endif
52
53   INTEGER    ::   nadv   ! choice of the formulation and scheme for the advection
54   LOGICAL    ::   lfirst=.TRUE.
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58
59CONTAINS
60
61   SUBROUTINE dyn_adv_tan( kt )
62      !!---------------------------------------------------------------------
63      !!                  ***  ROUTINE dyn_adv_tan  ***
64      !!
65      !! ** Purpose of the direct routine:
66      !!            compute the ocean momentum advection trend.
67      !!
68      !! ** Method  : - Update (ua,va) with the advection term following nadv
69      !!      NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T)
70      !!      a metric term is add to the coriolis term while in vector form
71      !!      it is the relative vorticity which is added to coriolis term
72      !!      (see dynvor module).
73      !!----------------------------------------------------------------------
74      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_tan')
78      !
79      SELECT CASE ( nadv )                     ! compute advection trend and add it to general trend
80      CASE ( 0 )
81                      CALL dyn_keg_tan     ( kt )    ! vector form : horizontal gradient of kinetic energy
82                      CALL dyn_zad_tan     ( kt )    ! vector form : vertical advection
83      CASE ( 1 )
84                      CALL dyn_adv_cen2_tan( kt )    ! 2nd order centered scheme
85      CASE ( 2 )
86                      CALL dyn_adv_ubs_tan ( kt )    ! 3rd order UBS      scheme
87      !
88      CASE (-1 )                                 ! esopa: test all possibility with control print
89                      CALL dyn_keg_tan     ( kt )
90                      CALL dyn_zad_tan     ( kt )
91                      CALL dyn_adv_cen2_tan( kt )
92                      CALL dyn_adv_ubs_tan ( kt )
93      END SELECT
94      !
95      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_tan')
96      !
97   END SUBROUTINE dyn_adv_tan
98
99   SUBROUTINE dyn_adv_adj( kt )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dyn_adv_adj  ***
102      !!
103      !! ** Purpose of the direct routine:
104      !!            compute the ocean momentum advection trend.
105      !!
106      !! ** Method  : - Update (ua,va) with the advection term following nadv
107      !!      NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T)
108      !!      a metric term is add to the coriolis term while in vector form
109      !!      it is the relative vorticity which is added to coriolis term
110      !!      (see dynvor module).
111      !!----------------------------------------------------------------------
112      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
113      !!----------------------------------------------------------------------
114      !
115      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_adj')
116      !
117      SELECT CASE ( nadv )                     ! compute advection trend and add it to general trend
118      CASE ( 0 )
119                      CALL dyn_zad_adj     ( kt )    ! vector form : vertical advection
120                      CALL dyn_keg_adj     ( kt )    ! vector form : horizontal gradient of kinetic energy
121      CASE ( 1 )
122         IF (lwp) WRITE(numout,*) 'dyn_adv_cen2_adj not available yet'
123         CALL abort
124      CASE ( 2 )
125         IF (lwp) WRITE(numout,*) 'dyn_adv_ubs_adj not available yet'
126         CALL abort
127      CASE (-1 )                                 ! esopa: test all possibility with control print
128                      CALL dyn_zad_adj     ( kt )
129                      CALL dyn_keg_adj     ( kt )
130      END SELECT
131      !
132      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_adj')
133      !
134   END SUBROUTINE dyn_adv_adj
135
136   SUBROUTINE dyn_adv_adj_tst( kumadt )
137      !!-----------------------------------------------------------------------
138      !!
139      !!                  ***  ROUTINE dyn_adv_adj_tst ***
140      !!
141      !! ** Purpose : Test the adjoint routine.
142      !!
143      !! ** Method  : Verify the scalar product
144      !!
145      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
146      !!
147      !!              where  L   = tangent routine
148      !!                     L^T = adjoint routine
149      !!                     W   = diagonal matrix of scale factors
150      !!                     dx  = input perturbation (random field)
151      !!                     dy  = L dx
152      !!
153      !!
154      !! History :
155      !!        ! 08-08 (A. Vidard)
156      !!-----------------------------------------------------------------------
157      !! * Modules used
158
159      !! * Arguments
160      INTEGER, INTENT(IN) :: &
161         & kumadt             ! Output unit
162
163      !! * Local declarations
164      INTEGER ::  &
165         & ji,    &        ! dummy loop indices
166         & jj,    &
167         & jk,    &
168         & jt
169
170      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
171         & zun_tlin,     & ! Tangent input:  now  u-velocity
172         & zvn_tlin,     & ! Tangent input:  now  v-velocity
173         & zwn_tlin,     & ! Tangent input:  now  w-velocity
174         & zua_tlin,     & ! Tangent input: after u-velocity
175         & zva_tlin,     & ! Tangent input: after u-velocity
176         & zua_tlout,    & ! Tangent output:after u-velocity
177         & zva_tlout,    & ! Tangent output:after v-velocity
178         & zua_adin,     & ! adjoint input: after u-velocity
179         & zva_adin,     & ! adjoint input: after v-velocity
180         & zun_adout,    & ! adjoint output: now  u-velocity
181         & zvn_adout,    & ! adjoint output: now  v-velocity
182         & zwn_adout,    & ! adjoint output: now  u-velocity
183         & zua_adout,    & ! adjoint output:after v-velocity
184         & zva_adout,    & ! adjoint output:after u-velocity
185         & zuvw            ! 3D random field for u, v and w
186
187      REAL(KIND=wp) ::   &
188         & zsp1,         & ! scalar product involving the tangent routine
189         & zsp1_1,       & !   scalar product components
190         & zsp1_2,       &
191         & zsp2,         & ! scalar product involving the adjoint routine
192         & zsp2_1,       & !   scalar product components
193         & zsp2_2,       &
194         & zsp2_3,       &
195         & zsp2_4,       &
196         & zsp2_5
197
198      CHARACTER(LEN=14) :: cl_name
199
200
201      ! Allocate memory
202
203      ALLOCATE( &
204         & zun_tlin(jpi,jpj,jpk),     &
205         & zvn_tlin(jpi,jpj,jpk),     &
206         & zwn_tlin(jpi,jpj,jpk),     &
207         & zua_tlin(jpi,jpj,jpk),     &
208         & zva_tlin(jpi,jpj,jpk),     &
209         & zua_tlout(jpi,jpj,jpk),    &
210         & zva_tlout(jpi,jpj,jpk),    &
211         & zua_adin(jpi,jpj,jpk),     &
212         & zva_adin(jpi,jpj,jpk),     &
213         & zun_adout(jpi,jpj,jpk),    &
214         & zvn_adout(jpi,jpj,jpk),    &
215         & zwn_adout(jpi,jpj,jpk),    &
216         & zua_adout(jpi,jpj,jpk),    &
217         & zva_adout(jpi,jpj,jpk),    &
218         & zuvw(jpi,jpj,jpk)           &
219         & )
220
221
222      !===================================================================================
223      ! 1)  dx = ( un_tl, vn_tl, ua_tl, va_tl )        --> dynkeg
224      ! and dy = ( ua_tl, va_tl )
225      ! 2)  dx = ( un_tl, vn_tl, wn_tl, ua_tl, va_tl ) --> dynkeg
226      ! and dy = ( ua_tl, va_tl )
227      ! 3)  dx = ( un_tl, vn_tl, wn_tl, ua_tl, va_tl ) --> dynadv
228      ! and dy = ( ua_tl, va_tl )
229      !======================================================================
230
231      DO jt = 1, 3
232      !--------------------------------------------------------------------
233      ! Reset the tangent and adjoint variables
234      !--------------------------------------------------------------------
235          zun_tlin(:,:,:)  = 0.0_wp
236          zvn_tlin(:,:,:)  = 0.0_wp
237          zwn_tlin(:,:,:)  = 0.0_wp
238          zua_tlin(:,:,:)  = 0.0_wp
239          zva_tlin(:,:,:)  = 0.0_wp
240          zua_tlout(:,:,:) = 0.0_wp
241          zva_tlout(:,:,:) = 0.0_wp
242          zua_adin(:,:,:)  = 0.0_wp
243          zva_adin(:,:,:)  = 0.0_wp
244          zun_adout(:,:,:) = 0.0_wp
245          zvn_adout(:,:,:) = 0.0_wp
246          zwn_adout(:,:,:) = 0.0_wp
247          zua_adout(:,:,:) = 0.0_wp
248          zva_adout(:,:,:) = 0.0_wp
249          zuvw(:,:,:)      = 0.0_wp
250
251          un_tl(:,:,:) = 0.0_wp
252          vn_tl(:,:,:) = 0.0_wp
253          wn_tl(:,:,:) = 0.0_wp
254          ua_tl(:,:,:) = 0.0_wp
255          va_tl(:,:,:) = 0.0_wp
256          un_ad(:,:,:) = 0.0_wp
257          vn_ad(:,:,:) = 0.0_wp
258          wn_ad(:,:,:) = 0.0_wp
259          ua_ad(:,:,:) = 0.0_wp
260          va_ad(:,:,:) = 0.0_wp
261
262      !--------------------------------------------------------------------
263      ! Initialize the tangent input with random noise: dx
264      !--------------------------------------------------------------------
265
266         CALL grid_random(  zuvw, 'U', 0.0_wp, stdu )
267         DO jk = 1, jpk
268            DO jj = nldj, nlej
269               DO ji = nldi, nlei
270                  zun_tlin(ji,jj,jk) = zuvw(ji,jj,jk)
271               END DO
272            END DO
273         END DO
274         CALL grid_random(  zuvw, 'V', 0.0_wp, stdv )
275         DO jk = 1, jpk
276            DO jj = nldj, nlej
277               DO ji = nldi, nlei
278                  zvn_tlin(ji,jj,jk) = zuvw(ji,jj,jk)
279               END DO
280            END DO
281         END DO
282
283         CALL grid_random(  zuvw, 'W', 0.0_wp, stdw )
284         DO jk = 1, jpk
285            DO jj = nldj, nlej
286               DO ji = nldi, nlei
287                  zwn_tlin(ji,jj,jk) = zuvw(ji,jj,jk)
288               END DO
289            END DO
290         END DO
291
292         CALL grid_random(  zuvw, 'U', 0.0_wp, stdu )
293         DO jk = 1, jpk
294            DO jj = nldj, nlej
295               DO ji = nldi, nlei
296                  zua_tlin(ji,jj,jk) = zuvw(ji,jj,jk)
297               END DO
298            END DO
299         END DO
300
301         CALL grid_random(  zuvw, 'V', 0.0_wp, stdv )
302         DO jk = 1, jpk
303            DO jj = nldj, nlej
304               DO ji = nldi, nlei
305                  zva_tlin(ji,jj,jk) = zuvw(ji,jj,jk)
306               END DO
307            END DO
308         END DO
309
310         un_tl(:,:,:) = zun_tlin(:,:,:)
311         vn_tl(:,:,:) = zvn_tlin(:,:,:)
312
313         wn_tl(:,:,:) = zwn_tlin(:,:,:)
314         ua_tl(:,:,:) = zua_tlin(:,:,:)
315         va_tl(:,:,:) = zva_tlin(:,:,:)
316
317         SELECT CASE ( jt )
318         CASE ( 1 )
319            CALL dyn_adv_init_tam
320            CALL dyn_keg_tan( nit000 )
321         CASE ( 2 )
322            CALL dyn_adv_init_tam
323            CALL dyn_zad_tan( nit000 )
324         CASE ( 3 )
325            CALL dyn_adv_tan ( nit000 )
326         END SELECT
327
328         zua_tlout(:,:,:) = ua_tl(:,:,:)
329         zva_tlout(:,:,:) = va_tl(:,:,:)
330
331         !--------------------------------------------------------------------
332         ! Initialize the adjoint variables: dy^* = W dy
333         !--------------------------------------------------------------------
334
335         DO jk = 1, jpk
336            DO jj = nldj, nlej
337               DO ji = nldi, nlei
338                  zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
339                       &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
340                       &               * umask(ji,jj,jk)
341                  zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
342                       &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
343                       &               * vmask(ji,jj,jk)
344               END DO
345            END DO
346         END DO
347         !--------------------------------------------------------------------
348         ! Compute the scalar product: ( L dx )^T W dy
349         !--------------------------------------------------------------------
350
351         zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin )
352         zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin )
353         zsp1   = zsp1_1 + zsp1_2
354
355         !--------------------------------------------------------------------
356         ! Call the adjoint routine: dx^* = L^T dy^*
357         !--------------------------------------------------------------------
358
359         ua_ad(:,:,:) = zua_adin(:,:,:)
360         va_ad(:,:,:) = zva_adin(:,:,:)
361
362         SELECT CASE ( jt )
363         CASE ( 1 )
364            CALL dyn_keg_adj( nitend )
365         CASE ( 2 )
366            CALL dyn_zad_adj( nitend )
367         CASE ( 3 )
368            CALL dyn_adv_adj( nitend )
369         END SELECT
370
371         zun_adout(:,:,:) = un_ad(:,:,:)
372         zvn_adout(:,:,:) = vn_ad(:,:,:)
373         zwn_adout(:,:,:) = wn_ad(:,:,:)
374         zua_adout(:,:,:) = ua_ad(:,:,:)
375         zva_adout(:,:,:) = va_ad(:,:,:)
376
377         zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
378         zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
379         IF (jt .EQ. 1) THEN
380            zsp2_3 = 0.0_wp
381         ELSE
382            zsp2_3 = DOT_PRODUCT( zwn_tlin, zwn_adout )
383         ENDIF
384         zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout )
385         zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout )
386         zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5
387
388         ! Compare the scalar products
389
390         SELECT CASE ( jt )
391         CASE ( 1 )
392            cl_name = 'dyn_keg_adj   '
393         CASE ( 2 )
394            cl_name = 'dyn_zad_adj   '
395         CASE ( 3 )
396            cl_name = 'dyn_adv_adj   '
397         END SELECT
398
399         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
400
401      ENDDO
402
403      DEALLOCATE( &
404         & zun_tlin,     &
405         & zvn_tlin,     &
406         & zwn_tlin,     &
407         & zua_tlin,     &
408         & zva_tlin,     &
409         & zua_tlout,    &
410         & zva_tlout,    &
411         & zua_adin,     &
412         & zva_adin,     &
413         & zun_adout,    &
414         & zvn_adout,    &
415         & zwn_adout,    &
416         & zua_adout,    &
417         & zva_adout,    &
418         & zuvw          &
419         & )
420
421   END SUBROUTINE dyn_adv_adj_tst
422  !!======================================================================
423   SUBROUTINE dyn_adv_init_tam
424      !!---------------------------------------------------------------------
425      !!                  ***  ROUTINE dyn_adv_ctl  ***
426      !!
427      !! ** Purpose :   Control the consistency between namelist options for
428      !!              momentum advection formulation & scheme and set nadv
429      !!----------------------------------------------------------------------
430      INTEGER ::   ioptio
431
432      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs
433      !!----------------------------------------------------------------------
434
435      IF (lfirst) THEN
436
437         REWIND ( numnam )               ! Read Namelist namdyn_adv : momentum advection scheme
438         READ   ( numnam, namdyn_adv )
439
440         IF(lwp) THEN                    ! Namelist print
441            WRITE(numout,*)
442            WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme'
443            WRITE(numout,*) '~~~~~~~~~~~'
444            WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
445            WRITE(numout,*) '          Vector/flux form (T/F)             ln_dynadv_vec  = ', ln_dynadv_vec
446            WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
447            WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
448         ENDIF
449
450         ioptio = 0                      ! Parameter control
451         IF( ln_dynadv_vec  )   ioptio = ioptio + 1
452         IF( ln_dynadv_cen2 )   ioptio = ioptio + 1
453         IF( ln_dynadv_ubs  )   ioptio = ioptio + 1
454         IF( lk_esopa       )   ioptio =          1
455
456         IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
457
458      !                               ! Set nadv
459         IF( ln_dynadv_vec  )   nadv =  0
460         IF( ln_dynadv_cen2 )   nadv =  1
461         IF( ln_dynadv_ubs  )   nadv =  2
462         IF( lk_esopa       )   nadv = -1
463
464         IF(lwp) THEN                    ! Print the choice
465            WRITE(numout,*)
466            IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
467            IF( nadv ==  1 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
468            IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
469            IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation'
470         ENDIF
471         !
472         lfirst = .FALSE.
473      END IF
474   END SUBROUTINE dyn_adv_init_tam
475#endif
476END MODULE dynadv_tam
Note: See TracBrowser for help on using the repository browser.