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 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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