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.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynzad.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.0 KB
Line 
1MODULE dynzad
2   !!======================================================================
3   !!                       ***  MODULE  dynzad  ***
4   !! Ocean dynamics : vertical advection trend
5   !!======================================================================
6   
7   !!----------------------------------------------------------------------
8   !!   dyn_zad      : vertical advection momentum trend
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers
12   USE dom_oce         ! ocean space and time domain
13   USE in_out_manager  ! I/O manager
14   USE trdmod          ! ocean dynamics trends
15   USE trdmod_oce      ! ocean variables trends
16   USE flxrnf          ! ocean runoffs
17
18   IMPLICIT NONE
19   PRIVATE
20   
21   !! * Accessibility
22   PUBLIC dyn_zad                ! routine called by step.F90
23
24   !! * Substitutions
25#  include "domzgr_substitute.h90"
26#  include "vectopt_loop_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LOCEAN-IPSL (2005)
29   !! $Header$
30   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
31   !!----------------------------------------------------------------------
32
33CONTAINS
34
35#if defined key_autotasking
36   !!----------------------------------------------------------------------
37   !!   'key_autotasking'                              j-k-i loops (j-slab)
38   !!----------------------------------------------------------------------
39
40   SUBROUTINE dyn_zad( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE dynzad  ***
43      !!
44      !! ** Purpose :   Compute the now vertical momentum advection trend and
45      !!      add it to the general trend of momentum equation.
46      !!
47      !! ** Method  :   Use j-slab (j-k-i loops) for auto-tasking
48      !!      The now vertical advection of momentum is given by:
49      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
50      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
51      !!      Add this trend to the general trend (ua,va):
52      !!         (ua,va) = (ua,va) + w dz(u,v)
53      !!
54      !! ** Action  : - Update (ua,va) with the vert. momentum advection trends
55      !!              - Save the trends in (utrd,vtrd) ('key_trddyn')
56      !!
57      !! History :
58      !!   6.0  !  91-01  (G. Madec) Original code
59      !!   7.0  !  91-11  (G. Madec)
60      !!   7.5  !  96-01  (G. Madec) statement function for e3
61      !!   8.5  !  02-07  (G. Madec) Free form, F90
62      !!   9.0  !  04-08  (C. Talandier) New trends organization
63      !!----------------------------------------------------------------------
64      !! * modules used
65      USE oce, ONLY:   zwuw => ta,   & ! use ta as 3D workspace
66                       zwvw => sa      ! use sa as 3D workspace
67
68      !! * Arguments
69      INTEGER, INTENT( in ) ::   kt    ! ocean time-step inedx
70     
71      !! * Local declarations
72      INTEGER  ::   ji, jj, jk         ! dummy loop indices
73      REAL(wp) ::   zvn, zua, zva      ! temporary scalars
74      REAL(wp), DIMENSION(jpi) ::   &
75         zww                           ! temporary workspace
76      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
77         ztdua, ztdva                  ! temporary workspace
78      !!----------------------------------------------------------------------
79     
80      IF( kt == nit000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
83         IF(lwp) WRITE(numout,*) '~~~~~~~   Auto-tasking case, j-slab, no vector opt.'
84      ENDIF
85
86      ! Save ua and va trends
87      IF( l_trddyn )   THEN
88         ztdua(:,:,:) = ua(:,:,:) 
89         ztdva(:,:,:) = va(:,:,:) 
90      ENDIF
91
92      !                                                ! ===============
93      DO jj = 2, jpjm1                                 !  Vertical slab
94         !                                             ! ===============
95
96         ! Vertical momentum advection at level w and u- and v- vertical
97         ! ----------------------------------------------------------------
98         DO jk = 2, jpkm1
99            ! vertical fluxes
100            DO ji = 2, jpi
101               zww(ji) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
102            END DO
103            ! vertical momentum advection at w-point
104            DO ji = 2, jpim1
105               zvn = 0.25 * e1t(ji,jj+1) * e2t(ji,jj+1) * wn(ji,jj+1,jk)
106               zwuw(ji,jj,jk) = ( zww(ji+1) + zww(ji) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
107               zwvw(ji,jj,jk) = ( zvn       + zww(ji) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
108            END DO 
109         END DO   
110
111         ! Surface and bottom values set to zero
112         DO ji = 2, jpim1
113            zwuw(ji,jj, 1 ) = 0.e0
114            zwvw(ji,jj, 1 ) = 0.e0
115            zwuw(ji,jj,jpk) = 0.e0
116            zwvw(ji,jj,jpk) = 0.e0
117         END DO 
118
119         ! Vertical momentum advection at u- and v-points
120         ! ----------------------------------------------
121         DO jk = 1, jpkm1
122            DO ji = 2, jpim1
123               ! vertical momentum advective trends
124               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
125               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
126               ! add the trends to the general momentum trends
127               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
128               va(ji,jj,jk) = va(ji,jj,jk) + zva
129            END DO 
130         END DO 
131         !                                             ! ===============
132      END DO                                           !   End of slab
133      !                                                ! ===============
134
135      ! save the vertical advection trends for diagnostic
136      ! momentum trends
137      IF( l_trddyn )   THEN
138         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
139         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
140
141         CALL trd_mod(ztdua, ztdva, jpdtdzad, 'DYN', kt)
142      ENDIF
143
144      IF(l_ctl) THEN         ! print sum trends (used for debugging)
145         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
146         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
147         WRITE(numout,*) ' zad  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
148         u_ctl = zua   ;   v_ctl = zva
149      ENDIF
150
151   END SUBROUTINE dyn_zad
152
153#else
154   !!----------------------------------------------------------------------
155   !!   Default option                             k-j-i loop (vector opt.)
156   !!----------------------------------------------------------------------
157
158   SUBROUTINE dyn_zad ( kt )
159      !!----------------------------------------------------------------------
160      !!                  ***  ROUTINE dynzad  ***
161      !!
162      !! ** Purpose :   Compute the now vertical momentum advection trend and
163      !!      add it to the general trend of momentum equation.
164      !!
165      !! ** Method  :   The now vertical advection of momentum is given by:
166      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
167      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
168      !!      Add this trend to the general trend (ua,va):
169      !!         (ua,va) = (ua,va) + w dz(u,v)
170      !!
171      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
172      !!              - Save the trends in (utrd,vtrd) ('key_trddyn')
173      !!
174      !! History :
175      !!   8.5  !  02-07  (G. Madec)  Original code
176      !!----------------------------------------------------------------------
177      !! * modules used
178      USE oce, ONLY:   zwuw => ta,   & ! use ta as 3D workspace
179                       zwvw => sa      ! use sa as 3D workspace
180      !! * Arguments
181      INTEGER, INTENT( in ) ::   kt    ! ocean time-step inedx
182     
183      !! * Local declarations
184      INTEGER  ::   ji, jj, jk         ! dummy loop indices
185      REAL(wp) ::   zua, zva           ! temporary scalars
186      REAL(wp), DIMENSION(jpi,jpj) ::   &
187         zww                           ! temporary  workspace
188      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
189         ztdua, ztdva                  ! temporary workspace
190      !!----------------------------------------------------------------------
191     
192      IF( kt == nit000 ) THEN
193         IF(lwp)WRITE(numout,*)
194         IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
195         IF(lwp)WRITE(numout,*) '~~~~~~~   vector optimization k-j-i loop'
196      ENDIF
197
198      ! Save ua and va trends
199      IF( l_trddyn )   THEN
200         ztdua(:,:,:) = ua(:,:,:) 
201         ztdva(:,:,:) = va(:,:,:) 
202      ENDIF
203     
204      ! Vertical momentum advection at level w and u- and v- vertical
205      ! -------------------------------------------------------------
206      DO jk = 2, jpkm1
207         ! vertical fluxes
208         DO jj = 2, jpj
209            DO ji = fs_2, jpi   ! vector opt.
210               zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
211            END DO
212         END DO
213         ! vertical momentum advection at w-point
214         DO jj = 2, jpjm1
215            DO ji = fs_2, fs_jpim1   ! vector opt.
216               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) )
217               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) )
218            END DO 
219         END DO   
220      END DO
221
222      ! Surface and bottom values set to zero
223      DO jj = 2, jpjm1
224         DO ji = fs_2, fs_jpim1   ! vector opt.
225            zwuw(ji,jj, 1 ) = 0.e0
226            zwvw(ji,jj, 1 ) = 0.e0
227            zwuw(ji,jj,jpk) = 0.e0
228            zwvw(ji,jj,jpk) = 0.e0
229         END DO 
230      END DO
231
232
233      ! Vertical momentum advection at u- and v-points
234      ! ----------------------------------------------
235      DO jk = 1, jpkm1
236         DO jj = 2, jpjm1
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238               ! vertical momentum advective trends
239               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
240               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
241               ! add the trends to the general momentum trends
242               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
243               va(ji,jj,jk) = va(ji,jj,jk) + zva
244            END DO 
245         END DO 
246      END DO
247
248      ! save the vertical advection trends for diagnostic
249      ! momentum trends
250      IF( l_trddyn )   THEN
251         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
252         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
253
254         CALL trd_mod(ztdua, ztdva, jpdtdzad, 'DYN', kt)
255      ENDIF
256
257      IF(l_ctl) THEN         ! print sum trends (used for debugging)
258         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
259         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
260         WRITE(numout,*) ' zad  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
261         u_ctl = zua   ;   v_ctl = zva
262      ENDIF
263
264   END SUBROUTINE dyn_zad
265#endif
266
267!!======================================================================
268END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.