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

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

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