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/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

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