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

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

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