source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90 @ 10789

Last change on this file since 10789 was 10789, checked in by davestorkey, 21 months ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: Convert first batch of DYN routines and "wn" → "ww".

  • Property svn:keywords set to Id
File size: 5.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   !!   NEMO     0.5  ! 2002-07  (G. Madec) Free form, F90
8   !!----------------------------------------------------------------------
9   
10   !!----------------------------------------------------------------------
11   !!   dyn_zad       : vertical advection momentum trend
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE sbc_oce        ! surface boundary condition: ocean
16   USE trd_oce        ! trends: ocean variables
17   USE trddyn         ! trend manager: dynamics
18   !
19   USE in_out_manager ! I/O manager
20   USE lib_mpp        ! MPP library
21   USE prtctl         ! Print control
22   USE timing         ! Timing
23
24   IMPLICIT NONE
25   PRIVATE
26   
27   PUBLIC   dyn_zad       ! routine called by dynadv.F90
28
29   !! * Substitutions
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE dyn_zad ( kt, ktlev, pu_rhs, pv_rhs )
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  :   The now vertical advection of momentum is given by:
46      !!         w dz(u) = pu_rhs + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(uu(:,:,:,ktlev)) ]
47      !!         w dz(v) = pv_rhs + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(vv(:,:,:,ktlev)) ]
48      !!      Add this trend to the general trend (pu_rhs,pv_rhs):
49      !!         (pu_rhs,pv_rhs) = (pu_rhs,pv_rhs) + w dz(u,v)
50      !!
51      !! ** Action  : - Update (pu_rhs,pv_rhs) with the vert. momentum adv. trends
52      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
53      !!----------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt             ! ocean time-step index
55      INTEGER, INTENT(in) ::   ktlev          ! time level index for source terms
56      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
57      !
58      INTEGER  ::   ji, jj, jk   ! dummy loop indices
59      REAL(wp) ::   zua, zva     ! local scalars
60      REAL(wp), DIMENSION(jpi,jpj)     ::   zww
61      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwuw, zwvw
62      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv
63      !!----------------------------------------------------------------------
64      !
65      IF( ln_timing )   CALL timing_start('dyn_zad')
66      !
67      IF( kt == nit000 ) THEN
68         IF(lwp) WRITE(numout,*)
69         IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
70      ENDIF
71
72      IF( l_trddyn )   THEN         ! Save pu_rhs and pv_rhs trends
73         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
74         ztrdu(:,:,:) = pu_rhs(:,:,:) 
75         ztrdv(:,:,:) = pv_rhs(:,:,:) 
76      ENDIF
77     
78      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
79         DO jj = 2, jpj                   ! vertical fluxes
80            DO ji = fs_2, jpi             ! vector opt.
81               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
82            END DO
83         END DO
84         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
85            DO ji = fs_2, fs_jpim1        ! vector opt.
86               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( uu(ji,jj,jk-1,ktlev) - uu(ji,jj,jk,ktlev) )
87               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vv(ji,jj,jk-1,ktlev) - vv(ji,jj,jk,ktlev) )
88            END DO 
89         END DO   
90      END DO
91      !
92      ! Surface and bottom advective fluxes set to zero
93      DO jj = 2, jpjm1       
94         DO ji = fs_2, fs_jpim1           ! vector opt.
95            zwuw(ji,jj, 1 ) = 0._wp
96            zwvw(ji,jj, 1 ) = 0._wp
97            zwuw(ji,jj,jpk) = 0._wp
98            zwvw(ji,jj,jpk) = 0._wp
99         END DO 
100      END DO
101      !
102      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points
103         DO jj = 2, jpjm1
104            DO ji = fs_2, fs_jpim1       ! vector opt.
105               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)
106               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)
107            END DO 
108         END DO 
109      END DO
110
111      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
112         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
113         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
114         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
115         DEALLOCATE( ztrdu, ztrdv ) 
116      ENDIF
117      !                             ! Control print
118      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
119         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
120      !
121      IF( ln_timing )   CALL timing_stop('dyn_zad')
122      !
123   END SUBROUTINE dyn_zad
124
125   !!======================================================================
126END MODULE dynzad
Note: See TracBrowser for help on using the repository browser.