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/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 8168

Last change on this file since 8168 was 8168, checked in by glong, 7 years ago

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