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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90 @ 4428

Last change on this file since 4428 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 8.9 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 :  2.0  ! 2006-08  (G. Madec, S. Theetten)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   dyn_adv_cen2       : flux form momentum advection (ln_dynadv_cen2=T)
13   !!                        trends using a 2nd order centred scheme 
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE trdmod_oce     ! ocean variables trends
18   USE trdmod         ! ocean dynamics trends
19   USE in_out_manager ! I/O manager
20   USE lib_mpp        ! MPP library
21   USE prtctl         ! Print control
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dyn_adv_cen2   ! routine called by step.F90
27
28   !! * Control permutation of array indices
29#  include "oce_ftrans.h90"
30#  include "dom_oce_ftrans.h90"
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dyn_adv_cen2( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dyn_adv_cen2  ***
45      !!
46      !! ** Purpose :   Compute the now momentum advection trend in flux form
47      !!              and the general trend of the momentum equation.
48      !!
49      !! ** Method  :   Trend evaluated using now fields (centered in time)
50      !!
51      !! ** Action  :   (ua,va) updated with the now vorticity term trend
52      !!----------------------------------------------------------------------
53      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
54      USE oce     , ONLY:   zfu   => ta       , zfv   => sa       ! (ta,sa) used as 3D workspace
55      USE wrk_nemo, ONLY:   zfu_t => wrk_3d_1 , zfv_t => wrk_3d_4 , zfu_uw =>wrk_3d_6   ! 3D workspaces
56      USE wrk_nemo, ONLY:   zfu_f => wrk_3d_2 , zfv_f => wrk_3d_5 , zfv_vw =>wrk_3d_7
57      USE wrk_nemo, ONLY:   zfw   => wrk_3d_3 
58      !! DCSE_NEMO: module variables renamed, need additional directives
59!FTRANS zfu :I :I :z
60!FTRANS zfv :I :I :z
61!FTRANS zfu_t :I :I :z
62!FTRANS zfv_t :I :I :z
63!FTRANS zfu_uw :I :I :z
64!FTRANS zfu_f :I :I :z
65!FTRANS zfv_f :I :I :z
66!FTRANS zfv_vw :I :I :z
67!FTRANS zfw :I :I :z
68      !
69      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
70      !
71      INTEGER  ::   ji, jj, jk   ! dummy loop indices
72      REAL(wp) ::   zbu, zbv     ! local scalars
73      !!----------------------------------------------------------------------
74
75      IF( kt == nit000 .AND. lwp ) THEN
76         WRITE(numout,*)
77         WRITE(numout,*) 'dyn_adv_cen2 : 2nd order flux form momentum advection'
78         WRITE(numout,*) '~~~~~~~~~~~~'
79      ENDIF
80
81      ! Check that global workspace arrays aren't already in use
82      IF( wrk_in_use(3, 1,2,3,4,5,6,7) ) THEN
83         CALL ctl_stop('dyn_adv_cen2 : requested workspace array unavailable')   ;   RETURN
84      ENDIF
85
86      IF( l_trddyn ) THEN           ! Save ua and va trends
87         zfu_uw(:,:,:) = ua(:,:,:)
88         zfv_vw(:,:,:) = va(:,:,:)
89      ENDIF
90
91      !                                      ! ====================== !
92      !                                      !  Horizontal advection  !
93      DO jk = 1, jpkm1                       ! ====================== !
94         !                                         ! horizontal volume fluxes
95         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
96         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
97         !
98         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
99            DO ji = 1, fs_jpim1   ! vector opt.
100               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
101               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) )
102               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) )
103               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
104            END DO
105         END DO
106         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
107            DO ji = fs_2, fs_jpim1   ! vector opt.
108               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
109               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
110               !
111               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
112                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
113               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
114                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
115            END DO
116         END DO
117      END DO
118      !
119      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
120         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
121         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
122         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
123         zfu_t(:,:,:) = ua(:,:,:)
124         zfv_t(:,:,:) = va(:,:,:)
125      ENDIF
126      !
127
128      !                                      ! ==================== !
129      !                                      !  Vertical advection  !
130      DO jk = 1, jpkm1                       ! ==================== !
131         !                                         ! Vertical volume fluxesÊ
132         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
133         !
134         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
135            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
136            zfv_vw(:,:,jpk) = 0.e0
137            !                                           ! Surface value :
138            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
139               zfu_uw(:,:, 1 ) = 0.e0   
140               zfv_vw(:,:, 1 ) = 0.e0
141            ELSE                                             ! constant volume : advection through the surface
142               DO jj = 2, jpjm1
143                  DO ji = fs_2, fs_jpim1
144                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
145                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
146                  END DO
147               END DO
148            ENDIF
149         ELSE                                      ! interior fluxes
150            DO jj = 2, jpjm1
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
153                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
154               END DO
155            END DO
156         ENDIF
157      END DO
158      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
159         DO jj = 2, jpjm1 
160            DO ji = fs_2, fs_jpim1   ! vector opt.
161               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
162                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
163               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
164                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
165            END DO
166         END DO
167      END DO
168      !
169      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
170         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
171         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
172         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
173      ENDIF
174      !                                            ! Control print
175      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask,   &
176         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
177      !
178      IF( wrk_not_released(3, 1,2,3,4,5,6,7) )   CALL ctl_stop('dyn_adv_cen2: failed to release workspace array')
179      !
180   END SUBROUTINE dyn_adv_cen2
181
182   !!==============================================================================
183END MODULE dynadv_cen2
Note: See TracBrowser for help on using the repository browser.