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

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN/dynzad_tam.F90 @ 1885

Last change on this file since 1885 was 1885, checked in by rblod, 14 years ago

add TAM sources
File size: 11.1 KB
Line 
1MODULE dynzad_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  dynzad_tam  ***
5   !! Ocean dynamics : vertical advection trend
6   !!                        Tangent and Adjoint module
7   !!======================================================================
8   !! History of the direct module:
9   !!            6.0  !  91-01  (G. Madec) Original code
10   !!            7.0  !  91-11  (G. Madec)
11   !!            7.5  !  96-01  (G. Madec) statement function for e3
12   !!            8.5  !  02-07  (G. Madec) j-k-i case: Original code
13   !!            8.5  !  02-07  (G. Madec) Free form, F90
14   !! History of the tam module:
15   !!            9.0  !  08-08  (A. Vidard) first version
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   dyn_zad_tan   : tangent of the vertical advection momentum trend
20   !!   dyn_zad_adj   : adjoint of the vertical advection momentum trend
21   !!----------------------------------------------------------------------
22   USE par_kind      , ONLY: & ! Precision variables
23      & wp
24   USE par_oce       , ONLY: & ! Ocean space and time domain variables
25      & jpi,                 &
26      & jpj,                 & 
27      & jpk,                 & 
28      & jpim1,               &
29      & jpjm1,               & 
30      & jpkm1
31   USE oce           , ONLY: & ! ocean dynamics and tracers
32      & un,                  &
33      & vn,                  &
34      & wn
35   USE dom_oce       , ONLY: & ! ocean space and time domain
36      & e1t,                 &
37      & e2t,                 &
38      & e1u,                 &
39      & e2u,                 &
40      & e1v,                 &
41      & e2v,                 &
42#if defined key_zco
43      & e3t_0
44#else
45      & e3u,                 &
46      & e3v
47#endif
48   USE in_out_manager, ONLY: & ! I/O manager
49      & lwp,                 &
50      & numout,              &
51      & nit000,              &
52      & nitend
53   USE oce_tam       , ONLY: & ! tam ocean dynamics and tracers
54      & un_tl,               &
55      & vn_tl,               &
56      & ua_tl,               &
57      & va_tl,               &
58      & un_ad,               &
59      & vn_ad,               &
60      & ua_ad,               &
61      & va_ad,               &
62      & wn_tl,               &
63      & wn_ad
64   IMPLICIT NONE
65   PRIVATE
66
67   PUBLIC   dyn_zad_tan   ! routine called by step_tam.F90
68   PUBLIC   dyn_zad_adj   ! routine called by step_tam.F90
69
70   !! * Substitutions
71#  include "domzgr_substitute.h90"
72#  include "vectopt_loop_substitute.h90"
73   !!----------------------------------------------------------------------
74   !!   OPA 9.0 , LOCEAN-IPSL (2005)
75   !! $Header$
76   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
77   !!----------------------------------------------------------------------
78
79CONTAINS
80   SUBROUTINE dyn_zad_tan ( kt )
81      !!----------------------------------------------------------------------
82      !!                  ***  ROUTINE dynzad_tan  ***
83      !!
84      !! ** Purpose of the direct routine:
85      !!      Compute the now vertical momentum advection trend and
86      !!      add it to the general trend of momentum equation.
87      !!
88      !! ** Method of the direct routine:
89      !!      The now vertical advection of momentum is given by:
90      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
91      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
92      !!      Add this trend to the general trend (ua,va):
93      !!         (ua,va) = (ua,va) + w dz(u,v)
94      !!
95      !! ** Action  : - Update (ua_tl,va_tl) with the vert. momentum advection trends
96      !!----------------------------------------------------------------------
97      !!
98      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
99      !!
100      INTEGER  ::   ji, jj, jk      ! dummy loop indices
101      REAL(wp) ::   zuatl, zvatl        ! temporary scalars
102      REAL(wp), DIMENSION(jpi,jpj)     ::   zww                ! 2D  workspace
103      REAL(wp), DIMENSION(jpi,jpj)     ::   zwwtl              ! 2D  workspace
104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuwtl, zwvwtl     ! 3D workspace
105      !!----------------------------------------------------------------------
106
107      IF( kt == nit000 ) THEN
108         IF(lwp)WRITE(numout,*)
109         IF(lwp)WRITE(numout,*) 'dyn_zad_tan : arakawa advection scheme'
110         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~'
111         CALL flush(numout)
112      ENDIF
113
114
115      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
116         DO jj = 2, jpj                   ! vertical fluxes
117            DO ji = fs_2, jpi             ! vector opt.
118               zww(  ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(   ji,jj,jk)
119               zwwtl(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn_tl(ji,jj,jk)
120            END DO
121         END DO
122         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
123            DO ji = fs_2, fs_jpim1        ! vector opt.
124               zwuwtl(ji,jj,jk) = ( zwwtl(ji+1,jj  ) + zwwtl(ji,jj) ) * ( un(   ji,jj,jk-1)-un(   ji,jj,jk) ) + &
125                  &               ( zww(  ji+1,jj  ) + zww(  ji,jj) ) * ( un_tl(ji,jj,jk-1)-un_tl(ji,jj,jk) )
126               zwvwtl(ji,jj,jk) = ( zwwtl(ji  ,jj+1) + zwwtl(ji,jj) ) * ( vn(   ji,jj,jk-1)-vn(   ji,jj,jk) ) + &
127                  &               ( zww(  ji  ,jj+1) + zww(  ji,jj) ) * ( vn_tl(ji,jj,jk-1)-vn_tl(ji,jj,jk) )
128            END DO
129         END DO
130      END DO
131      DO jj = 2, jpjm1              ! Surface and bottom values set to zero
132         DO ji = fs_2, fs_jpim1           ! vector opt.
133            zwuwtl(ji,jj, 1 ) = 0.e0
134            zwvwtl(ji,jj, 1 ) = 0.e0
135            zwuwtl(ji,jj,jpk) = 0.e0
136            zwvwtl(ji,jj,jpk) = 0.e0
137         END DO
138      END DO
139
140      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
141         DO jj = 2, jpjm1
142            DO ji = fs_2, fs_jpim1       ! vector opt.
143               !                         ! vertical momentum advective trends
144               zuatl = - ( zwuwtl(ji,jj,jk) + zwuwtl(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
145               zvatl = - ( zwvwtl(ji,jj,jk) + zwvwtl(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
146               !                         ! add the trends to the general momentum trends
147               ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zuatl
148               va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zvatl
149            END DO
150         END DO
151      END DO
152
153   END SUBROUTINE dyn_zad_tan
154
155   SUBROUTINE dyn_zad_adj ( kt )
156      !!----------------------------------------------------------------------
157      !!                  ***  ROUTINE dynzad_adj  ***
158      !!
159      !! ** Purpose of the direct routine:
160      !!      Compute the now vertical momentum advection trend and
161      !!      add it to the general trend of momentum equation.
162      !!
163      !! ** Method of the direct routine:
164      !!      The now vertical advection of momentum is given by:
165      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
166      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
167      !!      Add this trend to the general trend (ua,va):
168      !!         (ua,va) = (ua,va) + w dz(u,v)
169      !!
170      !! ** Action  : - Update (ua_tl,va_tl) with the vert. momentum advection trends
171      !!----------------------------------------------------------------------
172      !!
173      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
174      !!
175      INTEGER  ::   ji, jj, jk      ! dummy loop indices
176      REAL(wp) ::   zuaad, zvaad        ! temporary scalars
177      REAL(wp), DIMENSION(jpi,jpj)     ::   zww                ! 2D  workspace
178      REAL(wp), DIMENSION(jpi,jpj)     ::   zwwad              ! 2D  workspace
179      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuwad, zwvwad     ! 3D workspace
180      !!----------------------------------------------------------------------
181
182      IF( kt == nitend ) THEN
183         IF(lwp)WRITE(numout,*)
184         IF(lwp)WRITE(numout,*) 'dyn_zad_adj : arakawa advection scheme'
185         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~'
186      ENDIF
187
188      zuaad         = 0.0_wp
189      zvaad         = 0.0_wp
190      zwuwad(:,:,:) = 0.0_wp
191      zwvwad(:,:,:) = 0.0_wp
192      zwwad(:,:)    = 0.0_wp
193
194
195      DO jk = jpkm1, 1, -1             ! Vertical momentum advection at u- and v-points
196         DO jj = jpjm1, 2, -1
197            DO ji = fs_jpim1, fs_2, -1 ! vector opt.
198               !                       ! add the trends to the general momentum trends
199               zuaad = zuaad + ua_ad(ji,jj,jk)
200               zvaad = zvaad + va_ad(ji,jj,jk)
201               !                       ! vertical momentum advective trends
202               zwuwad(ji,jj,jk  ) = zwuwad(ji,jj,jk  ) - zuaad / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
203               zwuwad(ji,jj,jk+1) = zwuwad(ji,jj,jk+1) - zuaad / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
204               zuaad = 0.0_wp
205
206               zwvwad(ji,jj,jk  ) = zwvwad(ji,jj,jk  ) - zvaad / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
207               zwvwad(ji,jj,jk+1) = zwvwad(ji,jj,jk+1) - zvaad / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
208               zvaad = 0.0_wp
209            END DO
210         END DO
211      END DO
212      DO jk = jpkm1, 2, -1             ! Vertical momentum advection at level w and u- and v- vertical
213         DO jj = 2, jpj                   ! vertical fluxes
214            DO ji = fs_2, jpi             ! vector opt.
215               zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
216            END DO
217         END DO
218         DO jj = jpjm1, 2, -1          ! vertical momentum advection at w-point
219            DO ji = fs_jpim1, fs_2, -1 ! vector opt.
220               zwwad(ji,jj+1) = zwwad(ji,jj+1) + zwvwad(ji,jj,jk) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
221               zwwad(ji,jj  ) = zwwad(ji,jj  ) + zwvwad(ji,jj,jk) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
222               vn_ad(ji,jj,jk-1) = vn_ad(ji,jj,jk-1) + zwvwad(ji,jj,jk) * ( zww(ji,jj+1) + zww(ji,jj) )
223               vn_ad(ji,jj,jk  ) = vn_ad(ji,jj,jk  ) - zwvwad(ji,jj,jk) * ( zww(ji,jj+1) + zww(ji,jj) )
224               zwvwad(ji,jj,jk) = 0.0_wp
225
226               zwwad(ji+1,jj) = zwwad(ji+1,jj) + zwuwad(ji,jj,jk) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
227               zwwad(ji  ,jj) = zwwad(ji  ,jj) + zwuwad(ji,jj,jk) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
228               un_ad(ji,jj,jk-1) = un_ad(ji,jj,jk-1)  + zwuwad(ji,jj,jk) * ( zww(ji+1,jj) + zww(ji,jj) )
229               un_ad(ji,jj,jk  ) = un_ad(ji,jj,jk  )  - zwuwad(ji,jj,jk) * ( zww(ji+1,jj) + zww(ji,jj) )
230               zwuwad(ji,jj,jk) = 0.0_wp
231            END DO
232         END DO
233         DO jj = jpj, 2, -1                   ! vertical fluxes
234            DO ji = jpi, fs_2, -1             ! vector opt.
235               wn_ad(ji,jj,jk) = wn_ad(ji,jj,jk) + zwwad(ji,jj) * 0.25 * e1t(ji,jj) * e2t(ji,jj)
236               zwwad(ji,jj) = 0.0_wp
237            END DO
238         END DO
239      END DO
240
241   END SUBROUTINE dyn_zad_adj
242   SUBROUTINE dyn_zad_adj_tst( kumadt )         
243      INTEGER, INTENT(IN) :: &
244         & kumadt             ! Output unit
245      ! done in dynadv_tam
246   END SUBROUTINE dyn_zad_adj_tst
247#endif
248END MODULE dynzad_tam
Note: See TracBrowser for help on using the repository browser.