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/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA/traadv_tam.F90 @ 4583

Last change on this file since 4583 was 3611, checked in by pabouttier, 12 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007
  • Property svn:executable set to *
File size: 13.1 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   !!                 ! 12-07   (P.-A. Bouttier) Phase with 3.4 version
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   tra_adv      : compute ocean tracer advection trend
19   !!   tra_adv_init  : control the different options of advection scheme
20   !!----------------------------------------------------------------------
21   USE par_oce
22   USE oce_tam
23   USE oce
24   USE ldftra_oce
25   USE dom_oce         ! ocean space and time domain
26   USE traadv_cen2_tam
27   USE in_out_manager  ! I/O manager
28   USE prtctl          ! Print control
29   USE cla_tam
30   USE iom
31   USE lib_mpp
32   USE wrk_nemo
33   USE timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_adv_tan     ! routine called by steptan module
39   PUBLIC   tra_adv_adj     ! routine called by stepadj module
40   PUBLIC   tra_adv_init_tam ! routine called by stepadj module
41   PUBLIC   tra_adv_adj_tst ! routine called by tst module
42
43   !!* Namelist nam_traadv
44   LOGICAL, PUBLIC ::   ln_traadv_cen2   = .TRUE.       ! 2nd order centered scheme flag
45   LOGICAL, PUBLIC ::   ln_traadv_tvd    = .FALSE.      ! TVD scheme flag
46   LOGICAL, PUBLIC ::   ln_traadv_muscl  = .FALSE.      ! MUSCL scheme flag
47   LOGICAL, PUBLIC ::   ln_traadv_muscl2 = .FALSE.      ! MUSCL2 scheme flag
48   LOGICAL, PUBLIC ::   ln_traadv_ubs    = .FALSE.      ! UBS scheme flag
49   LOGICAL, PUBLIC ::   ln_traadv_qck    = .FALSE.      ! QUICKEST scheme flag
50
51   INTEGER ::   nadv   ! choice of the type of advection scheme
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56
57CONTAINS
58
59   SUBROUTINE tra_adv_tan( kt )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE tra_adv_tan  ***
62      !!
63      !! ** Purpose of the direct routine:
64      !!              compute the ocean tracer advection trend.
65      !!
66      !! ** Method  : - Update (ua,va) with the advection term following nadv
67      !!----------------------------------------------------------------------
68      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zuntl, zvntl, zwntl   ! effective velocity
69      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zun, zvn, zwn         ! effective velocity
70      !!
71      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
72      INTEGER               ::   jk   ! dummy loop index
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tan')
76      !
77      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
78      CALL wrk_alloc( jpi, jpj, jpk, zuntl, zvntl, zwntl )
79      !
80      IF( neuler == 0 .AND. kt == nit000 ) THEN     !at nit000
81         r2dtra(:) =  rdttra(:)                     !     = rdtra (restarting with Euler time stepping)
82      ELSEIF( kt <= nit000 + 1) THEN                !at nit000 or nit000+1
83         r2dtra(:) = 2._wp * rdttra(:)              !     = 2 rdttra (leapfrog)
84      ENDIF
85      !
86      IF( nn_cla == 1 )   CALL cla_traadv_tan( kt )
87      !
88      !                                               !==  effective transport  ==!
89      DO jk = 1, jpkm1
90         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
91         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
92         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
93         zuntl(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk)                  ! eulerian transport only
94         zvntl(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk)
95         zwntl(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn_tl(:,:,jk)
96      END DO
97      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
98      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
99      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
100      zuntl(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
101      zvntl(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
102      zwntl(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
103      !
104      IF(lwp) WRITE(numout,*) ' tra_adv_tam: 2nd order scheme is forced in TAM'
105      nadv = 1 ! force tra_adv_cen2 for tangent
106
107      SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend
108      CASE ( 1 ) ;
109         CALL tra_adv_cen2_tan ( kt, nit000, zun, zvn, zwn, tsn, zuntl, zvntl, zwntl, tsn_tl, tsa_tl, jpts )    ! 2nd order centered scheme
110      END SELECT
111      !
112      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tan')
113      !
114      CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn, zuntl, zvntl, zwntl )
115      !
116   END SUBROUTINE tra_adv_tan
117
118   SUBROUTINE tra_adv_adj( kt )
119      !!----------------------------------------------------------------------
120      !!                  ***  ROUTINE tra_adv_adj  ***
121      !!
122      !! ** Purpose of the direct routine:
123      !!              compute the ocean tracer advection trend.
124      !!
125      !! ** Method  : - Update (ua,va) with the advection term following nadv
126      !!----------------------------------------------------------------------
127      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zunad, zvnad, zwnad   ! effective velocity
128      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zun, zvn, zwn   ! effective velocity
129      INTEGER :: jk
130      !!
131      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
132      !!----------------------------------------------------------------------
133      IF( nn_timing == 1 )  CALL timing_start('tra_adv_adj')
134      !
135      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
136      CALL wrk_alloc( jpi, jpj, jpk, zunad, zvnad, zwnad )
137      !
138      zunad(:,:,:) = 0._wp                                                     ! no transport trough the bottom
139      zvnad(:,:,:) = 0._wp                                                     ! no transport trough the bottom
140      zwnad(:,:,:) = 0._wp
141      !
142      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
143         r2dtra(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
144      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
145         r2dtra(:) = 2._wp * rdttra(:)                   ! = 2 rdttra (leapfrog)
146      ENDIF
147      !
148      DO jk = 1, jpkm1
149         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
150         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
151         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
152      END DO
153      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
154      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
155      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
156      !
157      IF(lwp) WRITE(numout,*) ' tra_adv_tam: 2nd order scheme is forced in TAM'
158      nadv = 1 ! force tra_adv_cen2 for adjoint
159      !
160      SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend
161      CASE ( 1 ) ; CALL tra_adv_cen2_adj ( kt, nit000, zun, zvn, zwn, tsn, zunad, zvnad, zwnad, tsn_ad, tsa_ad, jpts )    ! 2nd order centered scheme
162      END SELECT
163      !
164      DO jk = jpkm1, 1, -1
165         un_ad(:,:,jk) = un_ad(:,:,jk) + e2u(:,:) * fse3u(:,:,jk) *  zunad(:,:,jk)
166         vn_ad(:,:,jk) = vn_ad(:,:,jk) + e1v(:,:) * fse3v(:,:,jk) *  zvnad(:,:,jk)
167         wn_ad(:,:,jk) = wn_ad(:,:,jk) + e1t(:,:) * e2t(:,:) *  zwnad(:,:,jk)
168      END DO
169      !
170      IF( nn_cla == 1 )   CALL cla_traadv_adj( kt )
171      !
172      !
173      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_adj')
174      !
175      CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn, zunad, zvnad, zwnad )
176      !
177   END SUBROUTINE tra_adv_adj
178   SUBROUTINE tra_adv_adj_tst( kumadt )
179      !!-----------------------------------------------------------------------
180      !!
181      !!                  ***  ROUTINE example_adj_tst ***
182      !!
183      !! ** Purpose : Test the adjoint routine.
184      !!
185      !! ** Method  : Verify the scalar product
186      !!
187      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
188      !!
189      !!              where  L   = tangent routine
190      !!                     L^T = adjoint routine
191      !!                     W   = diagonal matrix of scale factors
192      !!                     dx  = input perturbation (random field)
193      !!                     dy  = L dx
194      !!
195      !!
196      !! History :
197      !!        ! 08-08 (A. Vidard)
198      !!-----------------------------------------------------------------------
199      !! * Modules used
200
201      !! * Arguments
202      INTEGER, INTENT(IN) :: &
203         & kumadt             ! Output unit
204
205      CALL tra_adv_cen2_adj_tst( kumadt )
206
207   END SUBROUTINE tra_adv_adj_tst
208
209   SUBROUTINE tra_adv_init_tam
210      !!---------------------------------------------------------------------
211      !!                  ***  ROUTINE tra_adv_ctl_tam  ***
212      !!
213      !! ** Purpose :   Control the consistency between namelist options for
214      !!              tracer advection schemes and set nadv
215      !!----------------------------------------------------------------------
216      INTEGER ::   ioptio
217
218      NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd,    &
219         &                 ln_traadv_muscl, ln_traadv_muscl2, &
220         &                 ln_traadv_ubs  , ln_traadv_qck
221      !!----------------------------------------------------------------------
222
223!      REWIND ( numnam )               ! Read Namelist nam_traadv : tracer advection scheme
224!      READ   ( numnam, nam_traadv )
225
226      IF(lwp) THEN                    ! Namelist print
227         WRITE(numout,*)
228         WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme'
229         WRITE(numout,*) '~~~~~~~~~~~'
230         WRITE(numout,*) '       Namelist nam_traadv : chose a advection scheme for tracers'
231         WRITE(numout,*) '          2nd order advection scheme     ln_traadv_cen2   = ', ln_traadv_cen2
232         WRITE(numout,*) '          TVD advection scheme           ln_traadv_tvd    = ', ln_traadv_tvd
233         WRITE(numout,*) '          MUSCL  advection scheme        ln_traadv_muscl  = ', ln_traadv_muscl
234         WRITE(numout,*) '          MUSCL2 advection scheme        ln_traadv_muscl2 = ', ln_traadv_muscl2
235         WRITE(numout,*) '          UBS    advection scheme        ln_traadv_ubs    = ', ln_traadv_ubs
236         WRITE(numout,*) '          QUICKEST advection scheme      ln_traadv_qck    = ', ln_traadv_qck
237    ENDIF
238
239      ioptio = 0                      ! Parameter control
240      IF( ln_traadv_cen2   )   ioptio = ioptio + 1
241      IF( ln_traadv_tvd    )   ioptio = ioptio + 1
242      IF( ln_traadv_muscl  )   ioptio = ioptio + 1
243      IF( ln_traadv_muscl2 )   ioptio = ioptio + 1
244      IF( ln_traadv_ubs    )   ioptio = ioptio + 1
245      IF( ln_traadv_qck    )   ioptio = ioptio + 1
246      IF( lk_esopa         )   ioptio =          1
247
248      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist nam_traadv' )
249
250      IF( nn_cla == 1 .AND. .NOT. ln_traadv_cen2 )   &
251         &                CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' )
252
253      !                              ! Set nadv
254      IF( ln_traadv_cen2   )   nadv =  1
255      IF( ln_traadv_tvd    )   nadv =  2
256      IF( ln_traadv_muscl  )   nadv =  3
257      IF( ln_traadv_muscl2 )   nadv =  4
258      IF( ln_traadv_ubs    )   nadv =  5
259      IF( ln_traadv_qck    )   nadv =  6
260      IF( lk_esopa         )   nadv = -1
261
262      IF(lwp) THEN                   ! Print the choice
263         WRITE(numout,*)
264         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
265         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       Not Available in NEMO TAM'
266         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     Not Available in NEMO TAM'
267         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    Not Available in NEMO TAM'
268         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       Not Available in NEMO TAM'
269         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  Not Available in NEMO TAM'
270         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: Not Available in NEMO TAM'
271      ENDIF
272      !
273   END SUBROUTINE tra_adv_init_tam
274#endif
275
276
277  !!======================================================================
278END MODULE traadv_tam
Note: See TracBrowser for help on using the repository browser.