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

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

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

add TAM sources

File size: 12.1 KB
Line 
1MODULE dynldf_lap_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  dynldf_lap_tam  ***
5   !! Ocean dynamics:  lateral viscosity trend
6   !!                  Tangent and Adjoint Module
7   !!======================================================================
8
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_lap_tan  : update the momentum trend with the lateral diffusion
11   !!                      using an iso-level harmonic operator (tangent)
12   !!   dyn_ldf_lap_adj  : update the momentum trend with the lateral diffusion
13   !!                      using an iso-level harmonic operator (adjoint)
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE par_kind      , ONLY: & ! Precision variables
17      & wp
18   USE par_oce       , ONLY: & ! Ocean space and time domain variables
19      & jpi,                 &
20      & jpim1,               &
21      & jpjm1,               &
22      & jpkm1
23   USE oce_tam       , ONLY: &
24      & ua_tl,               &
25      & va_tl,               &
26      & ua_ad,               &
27      & va_ad,               &
28      & rotb_tl,             &
29      & hdivb_tl,            &
30      & rotb_ad,             &
31      & hdivb_ad
32   USE ldfdyn_oce    , ONLY: & ! ocean dynamics: lateral physics
33      & ahm2,                &
34      & ahm1,                &
35      & ahm0   
36   USE dom_oce       , ONLY: & ! Ocean space and time domain
37      & e1u,                 &
38      & e2u,                 &
39      & e1v,                 &
40      & e2v,                 &
41#if defined key_zco
42      & e3t_0
43#else
44      & e3u,                 &
45      & e3v,                 &
46      & e3f
47#endif
48   USE in_out_manager, ONLY: & ! I/O manager
49      & lwp,                 &
50      & numout,              & 
51      & nit000,              & 
52      & nitend
53
54   IMPLICIT NONE
55   PRIVATE
56
57   !! * Routine accessibility
58   PUBLIC dyn_ldf_lap_tan  ! called by dynldf_tam.F90
59   PUBLIC dyn_ldf_lap_adj  ! called by dynldf_tam.F90
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "ldfdyn_substitute.h90"
64#  include "vectopt_loop_substitute.h90"
65   !!----------------------------------------------------------------------
66
67CONTAINS
68
69   SUBROUTINE dyn_ldf_lap_tan( kt )
70      !!----------------------------------------------------------------------
71      !!                     ***  ROUTINE dyn_ldf_lap_tan  ***
72      !!                       
73      !! ** Purpose of the direct routine:
74      !!      Compute the before horizontal tracer (t & s) diffusive
75      !!      trend and add it to the general trend of tracer equation.
76      !!
77      !! ** Method of the direct routine:
78      !!      The before horizontal momentum diffusion trend is an
79      !!      harmonic operator (laplacian type) which separates the divergent
80      !!      and rotational parts of the flow.
81      !!      Its horizontal components are computed as follow:
82      !!         difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb]
83      !!         difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb]
84      !!      If lk_zco=T, e3f=e3u=e3v, the vertical scale factor are simplified
85      !!      in the rotational part of the diffusion.
86      !!      Add this before trend to the general trend (ua,va):
87      !!            (ua,va) = (ua,va) + (diffu,diffv)
88      !!      'key_trddyn' activated: the two components of the horizontal
89      !!                                 diffusion trend are saved.
90      !!
91      !! ** Action : - Update (ua,va) with the before iso-level harmonic
92      !!               mixing trend.
93      !!
94      !! History of the direct routine:
95      !!        !  90-09 (G. Madec) Original code
96      !!        !  91-11 (G. Madec)
97      !!        !  96-01 (G. Madec) statement function for e3 and ahm
98      !!   8.5  !  02-06 (G. Madec)  F90: Free form and module
99      !!   9.0  !  04-08 (C. Talandier) New trends organization
100      !! History of the tangent routine
101      !!   9.0  !  08-08 (A. Vidard) tangent of 9.0
102      !!----------------------------------------------------------------------
103      !! * Arguments
104      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
105      !! * Local declarations
106      INTEGER  ::   ji, jj, jk            ! dummy loop indices
107      REAL(wp) ::   &
108         zuatl, zvatl, ze2utl, ze1vtl             ! temporary scalars
109      !!----------------------------------------------------------------------
110
111      IF( kt == nit000 ) THEN
112         IF(lwp) WRITE(numout,*)
113         IF(lwp) WRITE(numout,*) 'dyn_ldf_tan: iso-level harmonic (laplacien) operator'
114         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
115      ENDIF
116
117      !                                                ! ===============
118      DO jk = 1, jpkm1                                 ! Horizontal slab
119         !                                             ! ===============
120         DO jj = 2, jpjm1
121            DO ji = fs_2, fs_jpim1   ! vector opt.
122# if defined key_zco
123               ! horizontal diffusive trends
124               ze2utl = rotb_tl (ji,jj,jk) * fsahmf(ji,jj,jk)
125               ze1vtl = hdivb_tl(ji,jj,jk) * fsahmt(ji,jj,jk)
126
127               zuatl = - ( ze2utl - rotb_tl (ji  ,jj-1,jk) * fsahmf(ji  ,jj-1,jk) ) / e2u(ji,jj)   &
128                  &    + ( hdivb_tl(ji+1,jj  ,jk) * fsahmt(ji+1,jj  ,jk) - ze1vtl ) / e1u(ji,jj)
129
130               zvatl = + ( ze2utl - rotb_tl (ji-1,jj  ,jk) * fsahmf(ji-1,jj  ,jk) ) / e1v(ji,jj)   &
131                  &    + ( hdivb_tl(ji  ,jj+1,jk) * fsahmt(ji  ,jj+1,jk) - ze1vtl ) / e2v(ji,jj)
132# else
133               ! horizontal diffusive trends
134               ze2utl = rotb_tl (ji,jj,jk) * fsahmf(ji,jj,jk) * fse3f(ji,jj,jk)
135               ze1vtl = hdivb_tl(ji,jj,jk) * fsahmt(ji,jj,jk)
136
137               zuatl = - ( ze2utl - rotb_tl(ji  ,jj-1,jk) * fsahmf(ji  ,jj-1,jk) * fse3f(ji  ,jj-1,jk) ) &
138                  &     / ( e2u(ji,jj) * fse3u(ji,jj,jk) )  &
139                  &    + ( hdivb_tl(ji+1,jj  ,jk) * fsahmt(ji+1,jj  ,jk) - ze1vtl ) / e1u(ji,jj)
140
141               zvatl = + ( ze2utl - rotb_tl(ji-1,jj  ,jk) * fsahmf(ji-1,jj  ,jk) * fse3f(ji-1,jj  ,jk) ) &
142                  &     / ( e1v(ji,jj) * fse3v(ji,jj,jk) )  &
143                  &    + ( hdivb_tl(ji  ,jj+1,jk) * fsahmt(ji  ,jj+1,jk) - ze1vtl ) / e2v(ji,jj)
144# endif
145               ! add it to the general momentum trends
146               ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zuatl
147               va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zvatl
148            END DO
149         END DO
150         !                                             ! ===============
151      END DO                                           !   End of slab
152      !                                                ! ===============
153
154   END SUBROUTINE dyn_ldf_lap_tan
155
156   SUBROUTINE dyn_ldf_lap_adj( kt )
157      !!----------------------------------------------------------------------
158      !!                     ***  ROUTINE dyn_ldf_lap_adj  ***
159      !!                       
160      !! ** Purpose of the direct routine:
161      !!      Compute the before horizontal tracer (t & s) diffusive
162      !!      trend and add it to the general trend of tracer equation.
163      !!
164      !! ** Method of the direct routine:
165      !!      The before horizontal momentum diffusion trend is an
166      !!      harmonic operator (laplacian type) which separates the divergent
167      !!      and rotational parts of the flow.
168      !!      Its horizontal components are computed as follow:
169      !!         difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb]
170      !!         difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb]
171      !!      If lk_zco=T, e3f=e3u=e3v, the vertical scale factor are simplified
172      !!      in the rotational part of the diffusion.
173      !!      Add this before trend to the general trend (ua,va):
174      !!            (ua,va) = (ua,va) + (diffu,diffv)
175      !!      'key_trddyn' activated: the two components of the horizontal
176      !!                                 diffusion trend are saved.
177      !!
178      !! ** Action : - Update (ua,va) with the before iso-level harmonic
179      !!               mixing trend.
180      !!
181      !! History of the direct routine:
182      !!        !  90-09 (G. Madec) Original code
183      !!        !  91-11 (G. Madec)
184      !!        !  96-01 (G. Madec) statement function for e3 and ahm
185      !!   8.5  !  02-06 (G. Madec)  F90: Free form and module
186      !!   9.0  !  04-08 (C. Talandier) New trends organization
187      !! History of the adjoint routine
188      !!   9.0  !  08-08 (A. Vidard) adjoint of 9.0
189      !!    -   !  09-01 (A. Weaver) misc. bug fixes and reorganization
190      !!----------------------------------------------------------------------
191      !! * Arguments
192      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
193      !! * Local declarations
194      INTEGER  ::   ji, jj, jk            ! dummy loop indices
195      REAL(wp) ::   &
196           zuaad , zvaad , ze2uad, ze1vad, &        ! temporary scalars
197         & zuaad1, zvaad1, zuaad2, zvaad2           ! temporary scalars
198      !!----------------------------------------------------------------------
199
200      IF( kt == nitend ) THEN
201         IF(lwp) WRITE(numout,*)
202         IF(lwp) WRITE(numout,*) 'dyn_ldf_adj: iso-level harmonic (laplacien) operator'
203         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
204      ENDIF
205
206      !                                                ! ===============
207      DO jk = jpkm1, 1, -1                                ! Horizontal slab
208         !                                             ! ===============
209         DO jj = jpjm1, 2, -1
210            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
211
212               ! add it to the general momentum trends
213               zuaad = ua_ad(ji,jj,jk) 
214               zvaad = va_ad(ji,jj,jk)
215# if defined key_zco
216               ! horizontal diffusive trends
217               zvaad1 = zvaad / e2v(ji,jj)
218               zvaad2 = zvaad / e1v(ji,jj)
219               zuaad1 = zuaad / e1u(ji,jj)
220               zuaad2 = zuaad / e2u(ji,jj)
221
222               ze1vad = - zvaad1 - zuaad1
223               ze2uad =   zvaad2 - zuaad2
224
225               rotb_ad (ji-1,jj  ,jk) = rotb_ad (ji-1,jj  ,jk) - zvaad2 * fsahmf(ji-1,jj  ,jk)
226               rotb_ad (ji  ,jj-1,jk) = rotb_ad (ji  ,jj-1,jk) + zuaad2 * fsahmf(ji  ,jj-1,jk)
227               rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) + ze2uad * fsahmf(ji  ,jj  ,jk)
228
229               hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) + zvaad1 * fsahmt(ji  ,jj+1,jk)
230               hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) + ze1vad * fsahmt(ji  ,jj  ,jk)
231               hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) + zuaad1 * fsahmt(ji+1,jj  ,jk)
232# else
233               ! horizontal diffusive trends
234               zvaad1 = zvaad /   e2v(ji,jj)
235               zvaad2 = zvaad / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) 
236               zuaad1 = zuaad /   e1u(ji,jj)
237               zuaad2 = zuaad / ( e2u(ji,jj) * fse3u(ji,jj,jk) )
238
239               ze1vad = - zvaad1 - zuaad1
240               ze2uad =   zvaad2 - zuaad2
241
242               rotb_ad (ji-1,jj  ,jk) = rotb_ad (ji-1,jj  ,jk) &
243                  &                   - zvaad2 * fsahmf(ji-1,jj  ,jk) * fse3f(ji-1,jj  ,jk)
244               rotb_ad (ji  ,jj-1,jk) = rotb_ad (ji  ,jj-1,jk) &
245                  &                   + zuaad2 * fsahmf(ji  ,jj-1,jk) * fse3f(ji  ,jj-1,jk)
246               rotb_ad (ji  ,jj  ,jk) = rotb_ad (ji  ,jj  ,jk) &
247                  &                   + ze2uad * fsahmf(ji  ,jj  ,jk) * fse3f(ji  ,jj  ,jk)
248
249               hdivb_ad(ji  ,jj+1,jk) = hdivb_ad(ji  ,jj+1,jk) &
250                  &                   + zvaad1 * fsahmt(ji  ,jj+1,jk)
251               hdivb_ad(ji  ,jj  ,jk) = hdivb_ad(ji  ,jj  ,jk) &
252                  &                   + ze1vad * fsahmt(ji  ,jj  ,jk)
253               hdivb_ad(ji+1,jj  ,jk) = hdivb_ad(ji+1,jj  ,jk) &
254                  &                   + zuaad1 * fsahmt(ji+1,jj  ,jk)
255
256# endif
257
258            END DO
259         END DO
260         !                                             ! ===============
261      END DO                                           !   End of slab
262      !                                                ! ===============
263
264   END SUBROUTINE dyn_ldf_lap_adj
265
266   !!======================================================================
267#endif
268END MODULE dynldf_lap_tam
Note: See TracBrowser for help on using the repository browser.