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.
traadv_tam.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA/traadv_tam.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago

refer to ticket #798

File size: 14.7 KB
Line 
1MODULE traadv_tam
2#if defined key_tam
3   !!==============================================================================
4   !!                       ***  MODULE  traadv_tam  ***
5   !! Ocean active tracers:  advection trend -
6   !!                        Tangent and Adjoint Module
7   !!==============================================================================
8   !! History of the direct routine:
9   !!            9.0  !  05-11  (G. Madec)  Original code
10   !! History of the TAM:
11   !!                 !  08-06  (A. Vidard) Skeleton
12   !!                 !  09-03  (F. Vigilant) Add tra_adv_ctl_tam routine
13   !!                                         Allow call to tra_eiv_tam/adj
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   tra_adv      : compute ocean tracer advection trend
18   !!   tra_adv_ctl  : control the different options of advection scheme
19   !!----------------------------------------------------------------------
20   USE par_kind         , ONLY: & ! Precision variables
21      & wp
22   USE par_oce          , ONLY: & ! Ocean space and time domain variables
23      & jpi,                    &
24      & jpj,                    & 
25      & jpk   
26   USE oce_tam          , ONLY: & ! ocean dynamics and active tracers
27      & un_tl,                  &
28      & vn_tl,                  &
29      & wn_tl,                  &
30      & un_ad,                  &
31      & vn_ad,                  &
32      & wn_ad
33   USE oce              , ONLY: & ! ocean dynamics and active tracers
34      & un,                     &
35      & vn,                     &
36      & wn
37   USE ldftra_oce       , ONLY: &
38      & lk_traldf_eiv
39   USE trabbl           , ONLY: &
40      & lk_trabbl_adv 
41   USE dom_oce         ! ocean space and time domain
42   USE traadv_cen2_tam, ONLY: & ! 2nd order centered scheme (tra_adv_cen2   routine)
43      & tra_adv_cen2_tan,     &
44      & tra_adv_cen2_adj,     &
45#if defined key_tst_tlm
46      & tra_adv_cen2_tlm_tst, &
47#endif
48      & tra_adv_cen2_adj_tst
49   USE traadv_eiv_tam, ONLY: & ! advection trend - eddy induced velocity (tra_adv_eiv   routine)
50      & tra_adv_eiv_tan,     &
51      & tra_adv_eiv_adj
52   USE in_out_manager  ! I/O manager
53   USE prtctl          ! Print control
54 
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC   tra_adv_tan     ! routine called by steptan module
59   PUBLIC   tra_adv_adj     ! routine called by stepadj module
60   PUBLIC   tra_adv_ctl_tam ! routine called by stepadj module
61   PUBLIC   tra_adv_adj_tst ! routine called by tst module
62#if defined key_tst_tlm
63   PUBLIC   tra_adv_tlm_tst ! routine called by tst module
64#endif
65   !!* Namelist nam_traadv
66   LOGICAL, PUBLIC ::   ln_traadv_cen2   = .TRUE.       ! 2nd order centered scheme flag
67   LOGICAL, PUBLIC ::   ln_traadv_tvd    = .FALSE.      ! TVD scheme flag
68   LOGICAL, PUBLIC ::   ln_traadv_muscl  = .FALSE.      ! MUSCL scheme flag
69   LOGICAL, PUBLIC ::   ln_traadv_muscl2 = .FALSE.      ! MUSCL2 scheme flag
70   LOGICAL, PUBLIC ::   ln_traadv_ubs    = .FALSE.      ! UBS scheme flag
71   LOGICAL, PUBLIC ::   ln_traadv_qck    = .FALSE.      ! QUICKEST scheme flag
72
73   INTEGER ::   nadv   ! choice of the type of advection scheme
74
75   !! * Substitutions
76#  include "domzgr_substitute.h90"
77#  include "vectopt_loop_substitute.h90"
78
79CONTAINS
80
81   SUBROUTINE tra_adv_tan( kt )
82      !!----------------------------------------------------------------------
83      !!                  ***  ROUTINE tra_adv_tan  ***
84      !!
85      !! ** Purpose of the direct routine:   
86      !!              compute the ocean tracer advection trend.
87      !!
88      !! ** Method  : - Update (ua,va) with the advection term following nadv
89      !!----------------------------------------------------------------------
90#if ( defined key_trabbl_adv || defined key_traldf_eiv )
91      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zuntl, zvntl, zwntl   ! effective velocity
92      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn         ! effective velocity
93#else
94      USE oce_tam, ONLY :                       zuntl => un_tl       ! the effective velocity is the
95      USE oce_tam, ONLY :                       zvntl => vn_tl       ! Eulerian velocity
96      USE oce_tam, ONLY :                       zwntl => wn_tl       !
97      USE oce, ONLY :                       zun => un       ! the effective velocity is the
98      USE oce, ONLY :                       zvn => vn       ! Eulerian velocity
99      USE oce, ONLY :                       zwn => wn       !
100#endif
101      !!
102      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
103      !!----------------------------------------------------------------------
104
105
106      IF( kt == nit000 )   CALL tra_adv_ctl_tam          ! initialisation & control of options
107
108#if defined key_trabbl_adv
109      zuntl(:,:,:) = un_tl(:,:,:) - u_bbl_tl(:,:,:)          ! add the bbl velocity
110      zvntl(:,:,:) = vn_tl(:,:,:) - v_bbl_tl(:,:,:)
111      zwntl(:,:,:) = wn_tl(:,:,:) + w_bbl_tl(:,:,:)
112      zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)          ! add the bbl velocity
113      zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
114      zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:)
115#endif
116      IF( lk_traldf_eiv ) THEN                       ! commpute and add the eiv velocity
117         IF( .NOT. lk_trabbl_adv ) THEN
118            zuntl(:,:,:) = un_tl(:,:,:)
119            zvntl(:,:,:) = vn_tl(:,:,:)
120            zwntl(:,:,:) = wn_tl(:,:,:)
121            zun(:,:,:) = un(:,:,:)
122            zvn(:,:,:) = vn(:,:,:)
123            zwn(:,:,:) = wn(:,:,:)
124         ENDIF
125         CALL tra_adv_eiv_tan( kt, zuntl, zvntl, zwntl ) 
126      ENDIF
127      IF(lwp) WRITE(numout,*) ' tra_adv_tam: 2nd order scheme is forced in TAM'
128      nadv = 1 ! force tra_adv_cen2 for tangent
129     
130      SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend
131      CASE ( 1 ) ; 
132         CALL tra_adv_cen2_tan ( kt, zun, zvn, zwn, zuntl, zvntl, zwntl )    ! 2nd order centered scheme
133      END SELECT
134      !
135   END SUBROUTINE tra_adv_tan
136
137   SUBROUTINE tra_adv_adj( kt )
138      !!----------------------------------------------------------------------
139      !!                  ***  ROUTINE tra_adv_adj  ***
140      !!
141      !! ** Purpose of the direct routine:   
142      !!              compute the ocean tracer advection trend.
143      !!
144      !! ** Method  : - Update (ua,va) with the advection term following nadv
145      !!----------------------------------------------------------------------
146#if ( defined key_trabbl_adv || defined key_traldf_eiv )
147      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zunad, zvnad, zwnad   ! effective velocity
148      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn   ! effective velocity
149#else
150      USE oce_tam, ONLY :                   zunad => un_ad       ! the effective velocity is the
151      USE oce_tam, ONLY :                   zvnad => vn_ad       ! Eulerian velocity
152      USE oce_tam, ONLY :                   zwnad => wn_ad       !
153      USE oce, ONLY :                       zun => un       ! the effective velocity is the
154      USE oce, ONLY :                       zvn => vn       ! Eulerian velocity
155      USE oce, ONLY :                       zwn => wn       !
156#endif
157      !!
158      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
159      !!----------------------------------------------------------------------
160
161      IF( kt == nitend )   CALL tra_adv_ctl_tam
162
163#if ( defined key_trabbl_adv || defined key_traldf_eiv )
164      zunad(:,:,:) = 0.0_wp 
165      zvnad(:,:,:) = 0.0_wp
166      zwnad(:,:,:) = 0.0_wp
167#endif
168#if defined key_trabbl_adv
169      zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)          ! add the bbl velocity
170      zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
171      zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:)
172#endif
173      IF( lk_traldf_eiv ) THEN                       ! commpute and add the eiv velocity
174         IF( .NOT. lk_trabbl_adv ) THEN
175            zun(:,:,:) = un(:,:,:)
176            zvn(:,:,:) = vn(:,:,:)
177            zwn(:,:,:) = wn(:,:,:)
178         ENDIF
179      ENDIF
180
181     IF(lwp) WRITE(numout,*) ' tra_adv_tam: 2nd order scheme is forced in TAM'
182      nadv = 1 ! force tra_adv_cen2 for adjoint
183     
184      SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend
185      CASE ( 1 ) ; CALL tra_adv_cen2_adj ( kt, zun, zvn, zwn, zunad, zvnad, zwnad )    ! 2nd order centered scheme
186      END SELECT
187      !
188      IF( lk_traldf_eiv ) THEN                       ! commpute and add the eiv velocity
189         CALL tra_adv_eiv_adj( kt, zunad, zvnad, zwnad ) 
190         IF( .NOT. lk_trabbl_adv ) THEN
191            un_ad(:,:,:) = un_ad(:,:,:) + zunad(:,:,:)
192            zunad(:,:,:) = 0.0_wp
193            vn_ad(:,:,:) = vn_ad(:,:,:) + zvnad(:,:,:)
194            zvnad(:,:,:) = 0.0_wp
195            wn_ad(:,:,:) = wn_ad(:,:,:) + zwnad(:,:,:)
196            zwnad(:,:,:) = 0.0_wp
197         ENDIF
198      ENDIF
199# if defined key_trabbl_adv
200      un_ad(:,:,:)    = un_ad(:,:,:)    + zunad(:,:,:)
201      u_bbl_ad(:,:,:) = u_bbl_ad(:,:,:) - zunad(:,:,:)
202      zunad(:,:,:) = 0.0_wp
203      vn_ad(:,:,:)    = vn_ad(:,:,:)    + zvnad(:,:,:)
204      v_bbl_ad(:,:,:) = v_bbl_ad(:,:,:) - zvnad(:,:,:)
205      zvnad(:,:,:) = 0.0_wp
206      wn_ad(:,:,:)    = wn_ad(:,:,:)    + zwnad(:,:,:)
207      w_bbl_ad(:,:,:) = w_bbl_ad(:,:,:) + zwnad(:,:,:)
208      zwnad(:,:,:) = 0.0_wp
209# endif
210   END SUBROUTINE tra_adv_adj
211SUBROUTINE tra_adv_adj_tst( kumadt )
212      !!-----------------------------------------------------------------------
213      !!
214      !!                  ***  ROUTINE example_adj_tst ***
215      !!
216      !! ** Purpose : Test the adjoint routine.
217      !!
218      !! ** Method  : Verify the scalar product
219      !!           
220      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
221      !!
222      !!              where  L   = tangent routine
223      !!                     L^T = adjoint routine
224      !!                     W   = diagonal matrix of scale factors
225      !!                     dx  = input perturbation (random field)
226      !!                     dy  = L dx
227      !!
228      !!                   
229      !! History :
230      !!        ! 08-08 (A. Vidard)
231      !!-----------------------------------------------------------------------
232      !! * Modules used
233
234      !! * Arguments
235      INTEGER, INTENT(IN) :: &
236         & kumadt             ! Output unit
237
238      CALL tra_adv_cen2_adj_tst( kumadt )
239
240   END SUBROUTINE tra_adv_adj_tst
241
242   SUBROUTINE tra_adv_ctl_tam
243      !!---------------------------------------------------------------------
244      !!                  ***  ROUTINE tra_adv_ctl_tam  ***
245      !!               
246      !! ** Purpose :   Control the consistency between namelist options for
247      !!              tracer advection schemes and set nadv
248      !!----------------------------------------------------------------------
249      INTEGER ::   ioptio
250
251      NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd,    &
252         &                 ln_traadv_muscl, ln_traadv_muscl2, &
253         &                 ln_traadv_ubs  , ln_traadv_qck
254      !!----------------------------------------------------------------------
255
256!      REWIND ( numnam )               ! Read Namelist nam_traadv : tracer advection scheme
257!      READ   ( numnam, nam_traadv )
258
259      IF(lwp) THEN                    ! Namelist print
260         WRITE(numout,*)
261         WRITE(numout,*) 'tra_adv_ctl : choice/control of the tracer advection scheme'
262         WRITE(numout,*) '~~~~~~~~~~~'
263         WRITE(numout,*) '       Namelist nam_traadv : chose a advection scheme for tracers'
264         WRITE(numout,*) '          2nd order advection scheme     ln_traadv_cen2   = ', ln_traadv_cen2
265         WRITE(numout,*) '          TVD advection scheme           ln_traadv_tvd    = ', ln_traadv_tvd
266         WRITE(numout,*) '          MUSCL  advection scheme        ln_traadv_muscl  = ', ln_traadv_muscl
267         WRITE(numout,*) '          MUSCL2 advection scheme        ln_traadv_muscl2 = ', ln_traadv_muscl2
268         WRITE(numout,*) '          UBS    advection scheme        ln_traadv_ubs    = ', ln_traadv_ubs
269         WRITE(numout,*) '          QUICKEST advection scheme      ln_traadv_qck    = ', ln_traadv_qck
270    ENDIF
271
272      ioptio = 0                      ! Parameter control
273      IF( ln_traadv_cen2   )   ioptio = ioptio + 1
274      IF( ln_traadv_tvd    )   ioptio = ioptio + 1
275      IF( ln_traadv_muscl  )   ioptio = ioptio + 1
276      IF( ln_traadv_muscl2 )   ioptio = ioptio + 1
277      IF( ln_traadv_ubs    )   ioptio = ioptio + 1
278      IF( ln_traadv_qck    )   ioptio = ioptio + 1
279      IF( lk_esopa         )   ioptio =          1
280
281      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist nam_traadv' )
282
283      IF( n_cla == 1 .AND. .NOT. ln_traadv_cen2 )   &
284         &                CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' )
285
286      !                              ! Set nadv
287      IF( ln_traadv_cen2   )   nadv =  1
288      IF( ln_traadv_tvd    )   nadv =  2
289      IF( ln_traadv_muscl  )   nadv =  3
290      IF( ln_traadv_muscl2 )   nadv =  4
291      IF( ln_traadv_ubs    )   nadv =  5
292      IF( ln_traadv_qck    )   nadv =  6
293      IF( lk_esopa         )   nadv = -1
294
295      IF(lwp) THEN                   ! Print the choice
296         WRITE(numout,*)
297         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
298         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       Not Available in NEMO TAM'
299         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     Not Available in NEMO TAM'
300         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    Not Available in NEMO TAM'
301         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       Not Available in NEMO TAM'
302         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  Not Available in NEMO TAM'
303         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: Not Available in NEMO TAM'
304      ENDIF
305      !
306   END SUBROUTINE tra_adv_ctl_tam
307#if defined key_tst_tlm
308SUBROUTINE tra_adv_tlm_tst( kumadt )
309      !!-----------------------------------------------------------------------
310      !!
311      !!                  ***  ROUTINE tra_adv_tlm_tst ***
312      !!
313      !! ** Purpose : Test the tangent routine.
314      !!
315      !! ** Method  : Verify the tangent with Taylor expansion
316      !!           
317      !!                 M(x+hdx) = M(x) + L(hdx) + O(h^2)
318      !!
319      !!              where  L   = tangent routine
320      !!                     M   = direct routine
321      !!                     dx  = input perturbation (random field)
322      !!                     h   = ration on perturbation
323      !!                   
324      !! History :
325      !!        ! 09-08 (A. Vigilant)
326      !!-----------------------------------------------------------------------
327      !! * Modules used
328
329      !! * Arguments
330      INTEGER, INTENT(IN) :: &
331         & kumadt             ! Output unit
332
333      CALL tra_adv_cen2_tlm_tst( kumadt )
334
335   END SUBROUTINE tra_adv_tlm_tst
336#endif
337#endif
338
339
340  !!======================================================================
341END MODULE traadv_tam
Note: See TracBrowser for help on using the repository browser.