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 branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago

changes after review
  • Property svn:keywords set to Id
File size: 13.7 KB
Line 
1MODULE dynzad
2   !!======================================================================
3   !!                       ***  MODULE  dynzad  ***
4   !! Ocean dynamics : vertical advection trend
5   !!======================================================================
6   !! History :  OPA  ! 1991-01  (G. Madec) Original code
7   !!            7.0  ! 1991-11  (G. Madec)
8   !!            7.5  ! 1996-01  (G. Madec) statement function for e3
9   !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
10   !!----------------------------------------------------------------------
11   
12   !!----------------------------------------------------------------------
13   !!   dyn_zad       : vertical advection momentum trend
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE sbc_oce        ! surface boundary condition: ocean
18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
21   USE in_out_manager ! I/O manager
22   USE lib_mpp        ! MPP library
23   USE prtctl         ! Print control
24   USE wrk_nemo       ! Memory Allocation
25   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29   
30   PUBLIC   dyn_zad       ! routine called by dynadv.F90
31   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
32
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dyn_zad ( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dynzad  ***
45      !!
46      !! ** Purpose :   Compute the now vertical momentum advection trend and
47      !!      add it to the general trend of momentum equation.
48      !!
49      !! ** Method  :   The now vertical advection of momentum is given by:
50      !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ]
51      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ]
52      !!      Add this trend to the general trend (ua,va):
53      !!         (ua,va) = (ua,va) + w dz(u,v)
54      !!
55      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
56      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
57      !!----------------------------------------------------------------------
58      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
59      !
60      INTEGER  ::   ji, jj, jk      ! dummy loop indices
61      REAL(wp) ::   zua, zva        ! temporary scalars
62      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw
63      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww
64      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
65      !!----------------------------------------------------------------------
66      !
67      IF( nn_timing == 1 )  CALL timing_start('dyn_zad')
68      !
69      CALL wrk_alloc( jpi,jpj, zww ) 
70      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 
71      !
72      IF( kt == nit000 ) THEN
73         IF(lwp)WRITE(numout,*)
74         IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
75      ENDIF
76
77      IF( l_trddyn )   THEN         ! Save ua and va trends
78         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
79!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
80         DO jk = 1, jpk
81            DO jj = 1, jpj
82               DO ji = 1, jpi
83                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
84                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
85               END DO
86            END DO
87         END DO
88      ENDIF
89     
90!$OMP PARALLEL
91      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
92!$OMP DO schedule(static) private(jj, ji)
93         DO jj = 2, jpj                   ! vertical fluxes
94            DO ji = fs_2, jpi             ! vector opt.
95               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
96            END DO
97         END DO
98!$OMP DO schedule(static) private(jj, ji)
99         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
100            DO ji = fs_2, fs_jpim1        ! vector opt.
101               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
102               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
103            END DO 
104         END DO   
105      END DO
106!$OMP END PARALLEL
107      !
108      ! Surface and bottom advective fluxes set to zero
109      IF ( ln_isfcav ) THEN
110!$OMP PARALLEL DO schedule(static) private(jj, ji)
111         DO jj = 2, jpjm1
112            DO ji = fs_2, fs_jpim1           ! vector opt.
113               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
114               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
115               zwuw(ji,jj,jpk) = 0._wp
116               zwvw(ji,jj,jpk) = 0._wp
117            END DO
118         END DO
119      ELSE
120!$OMP PARALLEL DO schedule(static) private(jj, ji)
121         DO jj = 2, jpjm1       
122            DO ji = fs_2, fs_jpim1           ! vector opt.
123               zwuw(ji,jj, 1 ) = 0._wp
124               zwvw(ji,jj, 1 ) = 0._wp
125               zwuw(ji,jj,jpk) = 0._wp
126               zwvw(ji,jj,jpk) = 0._wp
127            END DO 
128         END DO
129      END IF
130
131!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zua, zva)
132      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
133         DO jj = 2, jpjm1
134            DO ji = fs_2, fs_jpim1       ! vector opt.
135               !                         ! vertical momentum advective trends
136               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
137               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
138               !                         ! add the trends to the general momentum trends
139               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
140               va(ji,jj,jk) = va(ji,jj,jk) + zva
141            END DO 
142         END DO 
143      END DO
144
145      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
146!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
147         DO jk = 1, jpk
148            DO jj = 1, jpj
149               DO ji = 1, jpi
150                  ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
151                  ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
152               END DO
153            END DO
154         END DO
155         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
156         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
157      ENDIF
158      !                             ! Control print
159      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
160         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
161      !
162      CALL wrk_dealloc( jpi,jpj, zww ) 
163      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
164      !
165      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
166      !
167   END SUBROUTINE dyn_zad
168
169
170   SUBROUTINE dyn_zad_zts ( kt )
171      !!----------------------------------------------------------------------
172      !!                  ***  ROUTINE dynzad_zts  ***
173      !!
174      !! ** Purpose :   Compute the now vertical momentum advection trend and
175      !!      add it to the general trend of momentum equation. This version
176      !!      uses sub-timesteps for improved numerical stability with small
177      !!      vertical grid sizes. This is especially relevant when using
178      !!      embedded ice with thin surface boxes.
179      !!
180      !! ** Method  :   The now vertical advection of momentum is given by:
181      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
182      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
183      !!      Add this trend to the general trend (ua,va):
184      !!         (ua,va) = (ua,va) + w dz(u,v)
185      !!
186      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
187      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
188      !!----------------------------------------------------------------------
189      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
190      !
191      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
192      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
193      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
194      REAL(wp) ::   zua, zva        ! temporary scalars
195      REAL(wp) ::   zr_rdt          ! temporary scalar
196      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
197      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
198      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
199      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
200      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
201      !!----------------------------------------------------------------------
202      !
203      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
204      !
205      CALL wrk_alloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
206      CALL wrk_alloc( jpi,jpj,jpk,3,   zus , zvs ) 
207      !
208      IF( kt == nit000 ) THEN
209         IF(lwp)WRITE(numout,*)
210         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
211      ENDIF
212
213      IF( l_trddyn )   THEN         ! Save ua and va trends
214         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
215         ztrdu(:,:,:) = ua(:,:,:) 
216         ztrdv(:,:,:) = va(:,:,:) 
217      ENDIF
218     
219      IF( neuler == 0 .AND. kt == nit000 ) THEN
220          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
221      ELSE
222          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
223      ENDIF
224     
225      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
226         DO jj = 2, jpj                   
227            DO ji = fs_2, jpi             ! vector opt.
228               zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
229            END DO
230         END DO
231      END DO
232
233      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
234         DO ji = fs_2, fs_jpim1           ! vector opt.
235 !!gm missing ISF boundary condition
236            zwuw(ji,jj, 1 ) = 0._wp
237            zwvw(ji,jj, 1 ) = 0._wp
238            zwuw(ji,jj,jpk) = 0._wp
239            zwvw(ji,jj,jpk) = 0._wp
240         END DO 
241      END DO
242
243! Start with before values and use sub timestepping to reach after values
244
245      zus(:,:,:,1) = ub(:,:,:)
246      zvs(:,:,:,1) = vb(:,:,:)
247
248      DO jl = 1, jnzts                   ! Start of sub timestepping loop
249
250         IF( jl == 1 ) THEN              ! Euler forward to kick things off
251           jtb = 1   ;   jtn = 1   ;   jta = 2
252           zts = z2dtzts
253         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
254           jtb = 1   ;   jtn = 2   ;   jta = 3
255           zts = 2._wp * z2dtzts
256         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
257           jtb = MOD(jtb,3) + 1
258           jtn = MOD(jtn,3) + 1
259           jta = MOD(jta,3) + 1
260         ENDIF
261
262         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
263            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
264               DO ji = fs_2, fs_jpim1        ! vector opt.
265                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)
266                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)
267               END DO 
268            END DO   
269         END DO
270         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1       ! vector opt.
273                  !                         ! vertical momentum advective trends
274                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
275                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
276                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
277                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
278               END DO 
279            END DO 
280         END DO
281
282      END DO      ! End of sub timestepping loop
283
284      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
285      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1       ! vector opt.
288               !                         ! vertical momentum advective trends
289               !                         ! add the trends to the general momentum trends
290               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
291               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
292            END DO 
293         END DO 
294      END DO
295
296      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
297         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
298         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
299         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
300         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
301      ENDIF
302      !                             ! Control print
303      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
304         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
305      !
306      CALL wrk_dealloc( jpi,jpj,jpk,     zwuw, zwvw, zww ) 
307      CALL wrk_dealloc( jpi,jpj,jpk,3,   zus , zvs ) 
308      !
309      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
310      !
311   END SUBROUTINE dyn_zad_zts
312
313   !!======================================================================
314END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.