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.
dynadv_cen2.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynadv_cen2.F90 @ 1129

Last change on this file since 1129 was 1129, checked in by ctlod, 16 years ago

trunk: compilation error when key_trddyn is active related to the dynamics advection flux formulation, see ticket: #211
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 KB
Line 
1MODULE dynadv_cen2
2   !!======================================================================
3   !!                       ***  MODULE  dynadv  ***
4   !! Ocean dynamics: Update the momentum trend with the flux form advection
5   !!                 using a 2nd order centred scheme
6   !!======================================================================
7   !! History :  9.0  !  06-08  (G. Madec, S. Theetten)  Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_adv_cen2       : flux form momentum advection (ln_dynadv_cen2=T)
12   !!                        trends using a 2nd order centred scheme 
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE dynspg_oce     ! surface pressure gradient
17   USE in_out_manager ! I/O manager
18   USE dynspg_rl      ! surface pressure gradient
19   USE trdmod         ! ocean dynamics trends
20   USE trdmod_oce     ! ocean variables trends
21   USE prtctl         ! Print control
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC dyn_adv_cen2                        ! routine called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !!   OPA 9.0 , LODYC-IPSL  (2006)
34   !! $Header$
35   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE dyn_adv_cen2( kt )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE dyn_adv_cen2  ***
43      !!
44      !! ** Purpose :   Compute the now momentum advection trend in flux form
45      !!      and the general trend of the momentum equation.
46      !!
47      !! ** Method  :   Trend evaluated using now fields (centered in time)
48      !!
49      !! ** Action : - Update (ua,va) with the now vorticity term trend
50      !!----------------------------------------------------------------------
51      USE oce, ONLY:   zfu => ta,   & ! use ta as 3D workspace
52                       zfv => sa      ! use sa as 3D workspace
53
54      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
55
56      INTEGER  ::   ji, jj, jk              ! dummy loop indices
57      REAL(wp) ::   zua, zva, zbu, zbv      ! temporary scalars
58      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfu_t, zfu_f, zfu_uw   ! 3D workspace
59      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfv_t, zfv_f, zfv_vw   !  "      "
60      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfw                    !  "      "
61      !!----------------------------------------------------------------------
62
63      IF( kt == nit000 ) THEN
64         IF(lwp) WRITE(numout,*)
65         IF(lwp) WRITE(numout,*) 'dyn_adv_cen2 : 2nd order flux form momentum advection'
66         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
67      ENDIF
68
69      IF( l_trddyn ) THEN           ! Save ua and va trends
70         zfu_uw(:,:,:) = ua(:,:,:)
71         zfv_vw(:,:,:) = va(:,:,:)
72      ENDIF
73
74      ! I. Horizontal advection
75      ! -----------------------
76      !                                                ! ===============
77      DO jk = 1, jpkm1                                 ! Horizontal slab
78         !                                             ! ===============
79         ! horizontal volume fluxes
80         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
81         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
82
83         ! horizontal momentum fluxes at T- and F-point
84         DO jj = 1, jpjm1
85            DO ji = 1, fs_jpim1   ! vector opt.
86               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
87               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) )
88               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) )
89               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
90            END DO
91         END DO
92
93         ! divergence of horizontal momentum fluxes
94         DO jj = 2, jpjm1
95            DO ji = fs_2, fs_jpim1   ! vector opt.
96               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
97               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
98               ! horizontal advective trends
99               zua = - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)   &
100                  &     + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
101               zva = - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)   &
102                  &     + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
103               ! add it to the general tracer trends
104               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
105               va(ji,jj,jk) = va(ji,jj,jk) + zva
106            END DO
107         END DO
108         !                                             ! ===============
109      END DO                                           !   End of slab
110      !                                                ! ===============
111
112      IF( l_trddyn ) THEN      ! save the horizontal advection trend for diagnostic
113         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
114         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
115         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
116      ENDIF
117      !
118
119      ! II. Vertical advection
120      ! ----------------------
121
122      IF( l_trddyn ) THEN           ! Save ua and va trends
123         zfu_t(:,:,:) = ua(:,:,:)
124         zfv_t(:,:,:) = va(:,:,:)
125      ENDIF
126
127      ! Second order centered tracer flux at w-point
128      DO jk = 1, jpkm1
129         ! Vertical volume fluxes                   
130         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
131         ! Vertical advective fluxes                   
132         IF( jk == 1 ) THEN
133            zfu_uw(:,:,jpk) = 0.e0         ! Bottom value : flux set to zero
134            zfv_vw(:,:,jpk) = 0.e0
135            !                              ! Surface value
136            IF( lk_dynspg_rl ) THEN             ! rigid lid : flux set to zero
137               zfu_uw(:,:, 1 ) = 0.e0   
138               zfv_vw(:,:, 1 ) = 0.e0
139            ELSE                                ! free surface-constant volume
140               DO jj = 2, jpjm1
141                  DO ji = fs_2, fs_jpim1
142                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
143                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
144                  END DO
145               END DO
146            ENDIF
147         ELSE
148            !                              ! interior fluxes
149            DO jj = 2, jpjm1
150               DO ji = fs_2, fs_jpim1   ! vector opt.
151                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
152                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
153               END DO
154            END DO
155         ENDIF
156      END DO
157
158      ! momentum flux divergence at u-, v-points added to the general trend
159      DO jk = 1, jpkm1
160         DO jj = 2, jpjm1 
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zua = - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
163                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
164               zva = - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
165                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
166               ! add it to the general tracer trends
167               ua(ji,jj,jk) =  ua(ji,jj,jk) + zua
168               va(ji,jj,jk) =  va(ji,jj,jk) + zva
169            END DO
170         END DO
171      END DO
172
173      IF( l_trddyn ) THEN      ! save the vertical advection trend for diagnostic
174         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
175         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
176         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
177      ENDIF
178
179      !                             ! Control print
180      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask,   &
181         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
182      !
183   END SUBROUTINE dyn_adv_cen2
184
185   !!==============================================================================
186END MODULE dynadv_cen2
Note: See TracBrowser for help on using the repository browser.