source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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