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_cen2_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/traadv_cen2_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: 29.6 KB
Line 
1MODULE traadv_cen2_tam
2#if defined key_tam
3   !!======================================================================
4   !!                     ***  MODULE  traadv_cen2_tam  ***
5   !! Ocean active tracers:  horizontal & vertical advective trend
6   !!                         Tangent and Adjoint module
7   !!======================================================================
8   !! History :  8.2  ! 2001-08  (G. Madec, E. Durand)  trahad+trazad=traadv
9   !!            1.0  ! 2002-06  (G. Madec)  F90: Free form and module
10   !!            9.0  ! 2004-08  (C. Talandier) New trends organization
11   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
12   !!            2.0  ! 2006-04  (R. Benshila, G. Madec) Step reorganization
13   !!             -   ! 2006-07  (G. madec)  add ups_orca_set routine
14   !!            3.2  ! 2009-07  (G. Madec) add avmb, avtb in restart for cen2 advection
15   !! History of the T&A module
16   !!            9.0  ! 2008-12  (A. Vidard) original version
17   !!            3.2  ! 2010-04  (F. Vigilant)
18   !!            3.4  ! 2012-07  (P.-A. Bouttier)
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_adv_cen2 : update the tracer trend with the horizontal and
23   !!                  vertical advection trends using a seconder order
24   !!   ups_orca_set : allow mixed upstream/centered scheme in specific
25   !!                  area (set for orca 2 and 4 only)
26   !!----------------------------------------------------------------------
27   USE par_oce
28   USE oce
29   USE oce_tam
30   USE dom_oce
31   USE trc_oce
32   USE gridrandom
33   USE dotprodfld
34   USE in_out_manager
35   USE zdf_oce
36   USE tstool_tam
37   USE paresp
38   USE lib_mpp
39   USE wrk_nemo
40   use timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   tra_adv_cen2_tan    ! routine called by traadv_tam.F90
46   PUBLIC   tra_adv_cen2_adj    ! routine called by traadv_tam.F90
47   PUBLIC   tra_adv_cen2_adj_tst! routine called by tst.F90
48
49   !! * Substitutions
50#  include "domzgr_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
54   !! $Id: traadv_cen2.F90 1201 2008-09-24 13:24:21Z rblod $
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE tra_adv_cen2_tan( kt, kit000, pun, pvn, pwn, ptn, &
61      &                          pun_tl, pvn_tl, pwn_tl, ptn_tl, pta_tl, kjpt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE tra_adv_cen2_tan  ***
64      !!
65      !! ** Purpose of the direct routine:
66      !!      Compute the now trend due to the advection of tracers
67      !!      and add it to the general trend of passive tracer equations.
68      !!
69      !! ** Method  :   The advection is evaluated by a second order centered
70      !!      scheme using now fields (leap-frog scheme). In specific areas
71      !!      (vicinity of major river mouths, some straits, or where tn is
72      !!      approaching the freezing point) it is mixed with an upstream
73      !!      scheme for stability reasons.
74      !!         Part 0 : compute the upstream / centered flag
75      !!                  (3D array, zind, defined at T-point (0<zind<1))
76      !!         Part I : horizontal advection
77      !!       * centered flux:
78      !!               zcenu = e2u*e3u  un  mi(tn)
79      !!               zcenv = e1v*e3v  vn  mj(tn)
80      !!       * upstream flux:
81      !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0]
82      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
83      !!       * mixed upstream / centered horizontal advection scheme
84      !!               zcofi = max(zind(i+1), zind(i))
85      !!               zcofj = max(zind(j+1), zind(j))
86      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
87      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
88      !!       * horizontal advective trend (divergence of the fluxes)
89      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
90      !!       * Add this trend now to the general trend of tracer (ta,sa):
91      !!              (ta,sa) = (ta,sa) + ( zta , zsa )
92      !!       * trend diagnostic ('key_trdtra' defined): the trend is
93      !!      saved for diagnostics. The trends saved is expressed as
94      !!      Uh.gradh(T), i.e.
95      !!                     save trend = zta + tn divn
96      !!         In addition, the advective trend in the two horizontal direc-
97      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
98      !!      equal to (in s-coordinates, and similarly in z-coord.):
99      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] )
100      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  }
101      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so
102      !!      they vanish from the expression of the flux and divergence.
103      !!
104      !!         Part II : vertical advection
105      !!      For temperature (idem for salinity) the advective trend is com-
106      !!      puted as follows :
107      !!            zta = 1/e3t dk+1[ zwz ]
108      !!      where the vertical advective flux, zwz, is given by :
109      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
110      !!      with
111      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
112      !!        zcenu = centered flux = wn * mk(tn)
113      !!         The surface boundary condition is :
114      !!      variable volume (lk_vvl = T) : zero advective flux
115      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * tn(:,:,1)
116      !!         Add this trend now to the general trend of tracer (ta,sa):
117      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
118      !!         Trend diagnostic ('key_trdtra' defined): the trend is
119      !!      saved for diagnostics. The trends saved is expressed as :
120      !!             save trend =  w.gradz(T) = zta - tn divn.
121      !!
122      !! ** Action :  - update (ta,sa) with the now advective tracer trends
123      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
124      !!----------------------------------------------------------------------
125      !!
126      INTEGER , INTENT(in)                         ::   kt       ! ocean time-step index
127      INTEGER , INTENT(in   )                      ::   kit000   ! first time step index
128      INTEGER , INTENT(in   )                      ::   kjpt   ! first time step index
129      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun_tl   ! ocean velocity u-component
130      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn_tl   ! ocean velocity v-component
131      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn_tl   ! ocean velocity w-component
132      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun      ! ocean velocity u-component
133      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn      ! ocean velocity v-component
134      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn      ! ocean velocity w-component
135      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk,kjpt) ::   ptn     ! ocean velocity v-component
136      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) ::   ptn_tl, pta_tl     ! ocean velocity w-component
137      !!
138      INTEGER  ::   ji, jj, jk, jn                       ! dummy loop indices
139      REAL(wp) ::   zbtr, zhw, zhwtl,                 &  ! temporary scalars
140         &          ze3tr, zfui  , zfuitl  ,          &  !    "         "
141         &          zfvj  , zfvjtl, ztra                 !    "         "
142      REAL(wp) ::   zice                                 !    -         -
143      REAL(wp), POINTER, DIMENSION(:,:)     ::   ztfreez      ! 2D workspace
144      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwztl ! 3D workspace
145      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwxtl, zwytl ! 3D workspace
146      !!----------------------------------------------------------------------
147      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2_tan')
148      !
149      CALL wrk_alloc( jpi, jpj, ztfreez )
150      CALL wrk_alloc( jpi, jpj, jpk, zwztl, zwxtl, zwytl )
151      !
152      zwztl = 0._wp ; zwxtl = 0._wp ; zwytl = 0._wp
153      !
154      IF( kt == kit000 ) THEN
155         IF(lwp) WRITE(numout,*)
156         IF(lwp) WRITE(numout,*) 'tra_adv_cen2_tan : 2nd order centered advection scheme'
157         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
158         IF(lwp) WRITE(numout,*)
159         !
160      ENDIF
161      !
162      ! I. Horizontal advection
163      !    ====================
164      !
165      DO jn = 1, kjpt
166         DO jk = 1, jpkm1
167            !                        ! Second order centered tracer flux at u- and v-points
168            DO jj = 1, jpjm1
169               DO ji = 1, fs_jpim1   ! vector opt.
170                  ! volume fluxes * 1/2
171                  zfuitl = 0.5 * pun_tl(ji,jj,jk)
172                  zfvjtl = 0.5 * pvn_tl(ji,jj,jk)
173                  zfui   = 0.5 * pun(   ji,jj,jk)
174                  zfvj   = 0.5 * pvn(   ji,jj,jk)
175                  !
176                  ! centered scheme
177                  zwxtl(ji,jj,jk) = zfuitl * ( ptn(   ji,jj,jk, jn) + ptn(   ji+1,jj  ,jk, jn) ) &
178                     &            + zfui   * ( ptn_tl(ji,jj,jk, jn) + ptn_tl(ji+1,jj  ,jk, jn) )
179                  zwytl(ji,jj,jk) = zfvjtl * ( ptn(   ji,jj,jk, jn) + ptn(   ji  ,jj+1,jk, jn) ) &
180                     &            + zfvj   * ( ptn_tl(ji,jj,jk, jn) + ptn_tl(ji  ,jj+1,jk, jn) )
181               END DO
182            END DO
183         END DO
184
185         ! "zonal" mean advective heat and salt transport
186         ! ----------------------------------------------
187
188         ! II. Vertical advection
189         !     ==================
190         !
191         zwztl(:,:,jpk) = 0.0_wp                              ! Bottom value : flux set to zero
192
193         IF( lk_vvl ) THEN
194            zwztl(:,:, 1 ) = 0.e0                         ! volume variable
195         ELSE
196            zwztl(:,:, 1 ) = pwn(:,:,1) * ptn_tl(:,:,1,jn) &   ! linear free surface
197               &           + pwn_tl(:,:,1) * ptn(:,:,1,jn)
198         ENDIF
199         !
200         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
201            DO jj = 2, jpjm1
202               DO ji = fs_2, fs_jpim1   ! vector opt.
203                  ! velocity * 1/2
204                  zhwtl = 0.5 * pwn_tl(ji,jj,jk)
205                  zhw   = 0.5 * pwn(   ji,jj,jk)
206                  ! centered scheme
207                  zwztl(ji,jj,jk) = zhwtl * ( ptn(   ji,jj,jk,jn) + ptn(   ji,jj,jk-1,jn) ) &
208                     &            + zhw   * ( ptn_tl(ji,jj,jk,jn) + ptn_tl(ji,jj,jk-1,jn) )
209               END DO
210            END DO
211         END DO
212         !
213         DO jk = 1, jpkm1            ! divergence of Tracer flux added to the general trend
214            DO jj = 2, jpjm1
215               DO ji = fs_2, fs_jpim1   ! vector opt.
216                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
217                  ! vertical advective trends
218                  ztra = - zbtr * (  zwxtl(ji,jj,jk) - zwxtl(ji-1,jj  ,jk  )   &
219                  &                + zwytl(ji,jj,jk) - zwytl(ji  ,jj-1,jk  )   &
220                  &                + zwztl(ji,jj,jk) - zwztl(ji  ,jj  ,jk+1)  )
221                  ! advective trends added to the general tracer trends
222                  pta_tl(ji,jj,jk,jn) = pta_tl(ji,jj,jk,jn) + ztra
223               END DO
224            END DO
225         END DO
226      END DO
227      !
228      CALL wrk_dealloc( jpi, jpj, ztfreez )
229      CALL wrk_dealloc( jpi, jpj, jpk, zwztl, zwxtl, zwytl )
230      !
231      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2_tan')
232      !
233      !
234   END SUBROUTINE tra_adv_cen2_tan
235
236   SUBROUTINE tra_adv_cen2_adj( kt, kit000, pun, pvn, pwn, ptn, &
237      &                          pun_ad, pvn_ad, pwn_ad, ptn_ad, pta_ad, kjpt )
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE tra_adv_cen2_adj  ***
240      !!
241      !! ** Purpose of the direct routine:
242      !!      Compute the now trend due to the advection of tracers
243      !!      and add it to the general trend of passive tracer equations.
244      !!
245      !! ** Method  :   The advection is evaluated by a second order centered
246      !!      scheme using now fields (leap-frog scheme). In specific areas
247      !!      (vicinity of major river mouths, some straits, or where tn is
248      !!      approaching the freezing point) it is mixed with an upstream
249      !!      scheme for stability reasons.
250      !!         Part 0 : compute the upstream / centered flag
251      !!                  (3D array, zind, defined at T-point (0<zind<1))
252      !!         Part I : horizontal advection
253      !!       * centered flux:
254      !!               zcenu = e2u*e3u  un  mi(tn)
255      !!               zcenv = e1v*e3v  vn  mj(tn)
256      !!       * upstream flux:
257      !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0]
258      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
259      !!       * mixed upstream / centered horizontal advection scheme
260      !!               zcofi = max(zind(i+1), zind(i))
261      !!               zcofj = max(zind(j+1), zind(j))
262      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
263      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
264      !!       * horizontal advective trend (divergence of the fluxes)
265      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
266      !!       * Add this trend now to the general trend of tracer (ta,sa):
267      !!              (ta,sa) = (ta,sa) + ( zta , zsa )
268      !!       * trend diagnostic ('key_trdtra' defined): the trend is
269      !!      saved for diagnostics. The trends saved is expressed as
270      !!      Uh.gradh(T), i.e.
271      !!                     save trend = zta + tn divn
272      !!         In addition, the advective trend in the two horizontal direc-
273      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
274      !!      equal to (in s-coordinates, and similarly in z-coord.):
275      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] )
276      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  }
277      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so
278      !!      they vanish from the expression of the flux and divergence.
279      !!
280      !!         Part II : vertical advection
281      !!      For temperature (idem for salinity) the advective trend is com-
282      !!      puted as follows :
283      !!            zta = 1/e3t dk+1[ zwz ]
284      !!      where the vertical advective flux, zwz, is given by :
285      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
286      !!      with
287      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
288      !!        zcenu = centered flux = wn * mk(tn)
289      !!         The surface boundary condition is :
290      !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux
291      !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1)
292      !!         Add this trend now to the general trend of tracer (ta,sa):
293      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
294      !!         Trend diagnostic ('key_trdtra' defined): the trend is
295      !!      saved for diagnostics. The trends saved is expressed as :
296      !!             save trend =  w.gradz(T) = zta - tn divn.
297      !!
298      !! ** Action :  - update (ta,sa) with the now advective tracer trends
299      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
300      !!----------------------------------------------------------------------
301      !!
302      INTEGER , INTENT(in)                         ::   kt          ! ocean time-step index
303      INTEGER , INTENT(in)                         ::   kit000
304      INTEGER , INTENT(in)                         ::   kjpt
305      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun_ad   ! ocean velocity u-component
306      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn_ad   ! ocean velocity v-component
307      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn_ad   ! ocean velocity w-component
308      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun      ! ocean velocity u-component
309      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn      ! ocean velocity v-component
310      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn      ! ocean velocity w-component
311      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk,kjpt) ::   ptn
312      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) ::   ptn_ad, pta_ad
313      !!
314      INTEGER  ::   ji, jj, jk, jn                       ! dummy loop indices
315      REAL(wp) ::   zbtr, zhw, zhwad,                 &  ! temporary scalars
316         &          ze3tr, zfui  , zfuiad  ,          &  !    "         "
317         &          zfvj , zfvjad                        !    "         "
318      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwzad ! 3D workspace
319      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwxad, zwyad ! 3D workspace
320      !!----------------------------------------------------------------------
321      !
322      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2_adj')
323      !
324      CALL wrk_alloc( jpi, jpj, jpk, zwzad, zwyad, zwxad )
325      !
326      zhwad = 0.0_wp  ;  zfuiad = 0.0_wp  ;  zfvjad = 0.0_wp
327      zwxad(:,:,:) = 0.0_wp ; zwyad(:,:,:) = 0.0_wp ; zwzad(:,:,:) = 0.0_wp
328
329      IF( kt == nitend ) THEN
330         IF(lwp) WRITE(numout,*)
331         IF(lwp) WRITE(numout,*) 'tra_adv_cen2_adj : 2nd order centered advection scheme'
332         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
333         IF(lwp) WRITE(numout,*)
334         !
335      ENDIF
336      ! II. Vertical advection
337      !     ==================
338      !
339      DO jn = 1, kjpt
340         DO jk = jpkm1, 1, -1            ! divergence of Tracer flux added to the general trend
341            DO jj = jpjm1, 2, -1
342               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
343                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
344                  zwxad(ji,jj,jk) = zwxad(ji,jj,jk) - zbtr * pta_ad(ji,jj,jk,jn)
345                  zwyad(ji,jj,jk) = zwyad(ji,jj,jk) - zbtr * pta_ad(ji,jj,jk,jn)
346                  zwzad(ji,jj,jk) = zwzad(ji,jj,jk) - zbtr * pta_ad(ji,jj,jk,jn)
347                  zwxad(ji-1,jj,jk) = zwxad(ji-1,jj,jk) + zbtr * pta_ad(ji,jj,jk,jn)
348                  zwyad(ji,jj-1,jk) = zwyad(ji,jj-1,jk) + zbtr * pta_ad(ji,jj,jk,jn)
349                  zwzad(ji,jj,jk+1) = zwzad(ji,jj,jk+1) + zbtr * pta_ad(ji,jj,jk,jn)
350               END DO
351            END DO
352         END DO
353         !
354         DO jk = jpk, 2, -1              ! Second order centered tracer flux at w-point
355            DO jj = jpjm1, 2, -1
356               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
357                  zhw   = 0.5 * pwn(   ji,jj,jk)
358                  ! centered scheme
359                  zhwad = zwzad(ji,jj,jk) * ( ptn(   ji,jj,jk,jn) + ptn(   ji,jj,jk-1,jn) )
360                  ptn_ad(ji,jj,jk,jn) = ptn_ad(ji,jj,jk,jn) + zhw * zwzad(ji,jj,jk)
361                  ptn_ad(ji,jj,jk-1,jn) = ptn_ad(ji,jj,jk-1,jn) + zhw * zwzad(ji,jj,jk)
362                  zwzad(ji,jj,jk) = 0.0_wp
363                  pwn_ad(ji,jj,jk) = pwn_ad(ji,jj,jk) + 0.5 * zhwad
364                  zhwad = 0._wp
365              END DO
366            END DO
367         END DO
368
369         IF( lk_vvl ) THEN                                         ! Surface value : zero in variable volume
370            zwzad(:,:, 1 ) = 0.0_wp
371         ELSE                                                      !               : linear free surface case
372            pwn_ad(:,:,1) = pwn_ad(:,:,1) + zwzad(:,:, 1 ) * ptn(:,:,1,jn)
373            ptn_ad(:,:,1,jn)  = ptn_ad(:,:,1,jn)  + zwzad(:,:, 1 ) * pwn(:,:,1)
374            zwzad(:,:, 1 ) = 0.0_wp
375         ENDIF
376         !
377         zwzad(:,:,jpk) = 0.0_wp                                   ! Bottom value  : flux set to zero
378         ! "zonal" mean advective heat and salt transport
379         ! ----------------------------------------------
380         ! I. Horizontal advective fluxes
381         !    ====================
382         !
383         DO jk = 1, jpkm1
384            !                        ! Second order centered tracer flux at u- and v-points
385            DO jj = jpjm1, 1, -1
386               DO ji = fs_jpim1, 1, -1   ! vector opt.
387                  ! volume fluxes * 1/2
388                  zfui = 0.5 * pun(ji,jj,jk)
389                  zfvj = 0.5 * pvn(ji,jj,jk)
390                  ! centered scheme
391                  zfvjad = zwyad(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
392                  zfuiad = zwxad(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
393                  ptn_ad(ji  ,jj  ,jk,jn) = ptn_ad(ji  ,jj  ,jk,jn) + zwyad(ji,jj,jk) * zfvj
394                  ptn_ad(ji  ,jj+1,jk,jn) = ptn_ad(ji  ,jj+1,jk,jn) + zwyad(ji,jj,jk) * zfvj
395                  ptn_ad(ji  ,jj  ,jk,jn) = ptn_ad(ji  ,jj  ,jk,jn) + zwxad(ji,jj,jk) * zfui
396                  ptn_ad(ji+1,jj  ,jk,jn) = ptn_ad(ji+1,jj  ,jk,jn) + zwxad(ji,jj,jk) * zfui
397                  zwyad(ji  ,jj  ,jk) = 0.0_wp
398                  zwxad(ji  ,jj  ,jk) = 0.0_wp
399                  ! volume fluxes * 1/2
400                  pun_ad(ji,jj,jk) = pun_ad(ji,jj,jk) + 0.5 * zfuiad
401                  pvn_ad(ji,jj,jk) = pvn_ad(ji,jj,jk) + 0.5 * zfvjad
402                  zfuiad = 0._wp
403                  zfvjad = 0._wp
404                END DO
405            END DO
406         END DO
407      END DO
408      !
409      CALL wrk_dealloc( jpi, jpj, jpk, zwzad, zwyad, zwxad )
410      !
411      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2_adj')
412      !
413   END SUBROUTINE tra_adv_cen2_adj
414SUBROUTINE tra_adv_cen2_adj_tst( kumadt )
415      !!-----------------------------------------------------------------------
416      !!
417      !!                  ***  ROUTINE tra_adv_cen2_adj_tst ***
418      !!
419      !! ** Purpose : Test the adjoint routine.
420      !!
421      !! ** Method  : Verify the scalar product
422      !!
423      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
424      !!
425      !!              where  L   = tangent routine
426      !!                     L^T = adjoint routine
427      !!                     W   = diagonal matrix of scale factors
428      !!                     dx  = input perturbation (random field)
429      !!                     dy  = L dx
430      !!
431      !!
432      !! History :
433      !!        ! 08-08 (A. Vidard)
434      !!-----------------------------------------------------------------------
435      !! * Modules used
436
437      !! * Arguments
438      INTEGER, INTENT(IN) :: &
439         & kumadt             ! Output unit
440
441      !! * Local declarations
442      INTEGER ::  &
443         & ji,    &        ! dummy loop indices
444         & jj,    &
445         & jk
446      INTEGER, DIMENSION(jpi,jpj) :: &
447         & iseed_2d        ! 2D seed for the random number generator
448      REAL(KIND=wp) :: &
449         & zsp1,         & ! scalar product involving the tangent routine
450         & zsp2            ! scalar product involving the adjoint routine
451      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
452         & zun_tlin ,     zvn_tlin ,     zwn_tlin ,     ztn_tlin ,     & ! Tangent input
453         & zsn_tlin ,     zta_tlin ,     zsa_tlin ,     & ! Tangent input
454         & zun_adout,     zvn_adout,     zwn_adout,     ztn_adout,     & ! Adjoint output
455         & zsn_adout,     zta_adout,     zsa_adout,     & ! Adjoint output
456         & zta_tlout,     zsa_tlout,     & ! Tangent output
457         & zta_adin ,     zsa_adin ,     & ! Adjoint input
458         & zr             ! 3D random field
459      CHARACTER(LEN=14) ::&
460         & cl_name
461      ! Allocate memory
462
463      ALLOCATE( &
464         & zun_tlin( jpi,jpj,jpk),     zvn_tlin( jpi,jpj,jpk),     zwn_tlin( jpi,jpj,jpk),     &
465         & ztn_tlin( jpi,jpj,jpk),     zsn_tlin( jpi,jpj,jpk),     zta_tlin( jpi,jpj,jpk),     &
466         & zsa_tlin( jpi,jpj,jpk),     zta_tlout(jpi,jpj,jpk),     zsa_tlout(jpi,jpj,jpk),     &
467         & zta_adin( jpi,jpj,jpk),     zsa_adin( jpi,jpj,jpk),     zun_adout(jpi,jpj,jpk),     &
468         & zvn_adout(jpi,jpj,jpk),     zwn_adout(jpi,jpj,jpk),     ztn_adout(jpi,jpj,jpk),     &
469         & zsn_adout(jpi,jpj,jpk),     zta_adout(jpi,jpj,jpk),     zsa_adout(jpi,jpj,jpk),     &
470         & zr(       jpi,jpj,jpk)      &
471         & )
472      !==================================================================
473      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
474      !    dy = ( hdivb_tl, hdivn_tl )
475      !==================================================================
476
477      !--------------------------------------------------------------------
478      ! Reset the tangent and adjoint variables
479      !--------------------------------------------------------------------
480      zun_tlin( :,:,:) = 0.0_wp
481      zvn_tlin( :,:,:) = 0.0_wp
482      zwn_tlin( :,:,:) = 0.0_wp
483      ztn_tlin( :,:,:) = 0.0_wp
484      zsn_tlin( :,:,:) = 0.0_wp
485      zta_tlin( :,:,:) = 0.0_wp
486      zsa_tlin( :,:,:) = 0.0_wp
487      zta_tlout(:,:,:) = 0.0_wp
488      zsa_tlout(:,:,:) = 0.0_wp
489      zta_adin( :,:,:) = 0.0_wp
490      zsa_adin( :,:,:) = 0.0_wp
491      zun_adout(:,:,:) = 0.0_wp
492      zvn_adout(:,:,:) = 0.0_wp
493      zwn_adout(:,:,:) = 0.0_wp
494      ztn_adout(:,:,:) = 0.0_wp
495      zsn_adout(:,:,:) = 0.0_wp
496      zta_adout(:,:,:) = 0.0_wp
497      zsa_adout(:,:,:) = 0.0_wp
498      zr(       :,:,:) = 0.0_wp
499
500      tsn_ad(:,:,:,jp_tem) = 0.0_wp
501      tsn_ad(:,:,:,jp_sal) = 0.0_wp
502
503
504      !--------------------------------------------------------------------
505      ! Initialize the tangent input with random noise: dx
506      !--------------------------------------------------------------------
507
508      CALL grid_random(  zr, 'U', 0.0_wp, stdu )
509      DO jk = 1, jpk
510        DO jj = nldj, nlej
511           DO ji = nldi, nlei
512              zun_tlin(ji,jj,jk) = zr(ji,jj,jk)
513            END DO
514         END DO
515      END DO
516      CALL grid_random(  zr, 'V', 0.0_wp, stdv )
517      DO jk = 1, jpk
518        DO jj = nldj, nlej
519           DO ji = nldi, nlei
520              zvn_tlin(ji,jj,jk) = zr(ji,jj,jk)
521            END DO
522         END DO
523      END DO
524      CALL grid_random(  zr, 'W', 0.0_wp, stdw )
525      DO jk = 1, jpk
526        DO jj = nldj, nlej
527           DO ji = nldi, nlei
528              zwn_tlin(ji,jj,jk) = zr(ji,jj,jk)
529            END DO
530         END DO
531      END DO
532      CALL grid_random(  zr, 'T', 0.0_wp, stdt )
533      DO jk = 1, jpk
534        DO jj = nldj, nlej
535           DO ji = nldi, nlei
536              ztn_tlin(ji,jj,jk) = zr(ji,jj,jk)
537            END DO
538         END DO
539      END DO
540      CALL grid_random(  zr, 'T', 0.0_wp, stds )
541      DO jk = 1, jpk
542        DO jj = nldj, nlej
543           DO ji = nldi, nlei
544              zsn_tlin(ji,jj,jk) = zr(ji,jj,jk)
545            END DO
546         END DO
547      END DO
548      CALL grid_random(  zr, 'T', 0.0_wp, stdt )
549      DO jk = 1, jpk
550        DO jj = nldj, nlej
551           DO ji = nldi, nlei
552              zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
553            END DO
554         END DO
555      END DO
556      CALL grid_random(  zr, 'T', 0.0_wp, stds )
557      DO jk = 1, jpk
558        DO jj = nldj, nlej
559           DO ji = nldi, nlei
560              zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
561            END DO
562         END DO
563      END DO
564
565      tsn_tl(:,:,:,jp_tem) = ztn_tlin(:,:,:)
566      tsn_tl(:,:,:,jp_sal) = zsn_tlin(:,:,:)
567      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
568      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
569
570      CALL tra_adv_cen2_tan(nit000, nit000, un, vn, wn, tsn, zun_tlin, zvn_tlin, zwn_tlin, tsn_tl, tsa_tl, 2)
571
572      zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
573      zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
574
575      !--------------------------------------------------------------------
576      ! Initialize the adjoint variables: dy^* = W dy
577      !--------------------------------------------------------------------
578
579      DO jk = 1, jpk
580        DO jj = nldj, nlej
581           DO ji = nldi, nlei
582              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
583                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
584                 &               * tmask(ji,jj,jk) * wesp_t(jk)
585              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
586                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
587                 &               * tmask(ji,jj,jk) * wesp_s(jk)
588            END DO
589         END DO
590      END DO
591      !--------------------------------------------------------------------
592      ! Compute the scalar product: ( L dx )^T W dy
593      !--------------------------------------------------------------------
594
595      zsp1 = DOT_PRODUCT( zsa_tlout, zsa_adin ) &
596         & + DOT_PRODUCT( zta_tlout, zta_adin )
597
598      !--------------------------------------------------------------------
599      ! Call the adjoint routine: dx^* = L^T dy^*
600      !--------------------------------------------------------------------
601      tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
602      tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:)
603
604      CALL tra_adv_cen2_adj(nit000, nit000, un, vn, wn, tsn, zun_adout, zvn_adout, zwn_adout, tsn_ad, tsa_ad, 2)
605
606      ztn_adout(:,:,:) = tsn_ad(:,:,:,jp_tem)
607      zsn_adout(:,:,:) = tsn_ad(:,:,:,jp_sal)
608      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
609      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
610
611      zsp2 = DOT_PRODUCT( zun_tlin, zun_adout ) &
612         & + DOT_PRODUCT( zvn_tlin, zvn_adout ) &
613         & + DOT_PRODUCT( zwn_tlin, zwn_adout ) &
614         & + DOT_PRODUCT( ztn_tlin, ztn_adout ) &
615         & + DOT_PRODUCT( zsn_tlin, zsn_adout ) &
616         & + DOT_PRODUCT( zta_tlin, zta_adout ) &
617         & + DOT_PRODUCT( zsa_tlin, zsa_adout )
618
619      ! 14 char:'12345678901234'
620      cl_name = 'tra_adv_cen2  '
621      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
622
623      DEALLOCATE(         &
624         & zun_tlin ,     zvn_tlin ,     zwn_tlin ,     ztn_tlin ,     zsn_tlin ,     &
625         & zta_tlin ,     zsa_tlin ,     zta_tlout,     zsa_tlout,     zta_adin ,     &
626         & zsa_adin ,     zun_adout,     zvn_adout,     zwn_adout,     ztn_adout,     &
627         & zsn_adout,     zta_adout,     zsa_adout,     zr                            &
628         & )
629
630   END SUBROUTINE tra_adv_cen2_adj_tst
631#endif
632   !!======================================================================
633END MODULE traadv_cen2_tam
Note: See TracBrowser for help on using the repository browser.