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 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 13 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
RevLine 
[643]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   !!======================================================================
[1566]7   !! History :  2.0  ! 2006-08  (G. Madec, S. Theetten)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
[643]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
[1566]17   USE trdmod_oce     ! ocean variables trends
18   USE trdmod         ! ocean dynamics trends
[643]19   USE in_out_manager ! I/O manager
[2715]20   USE lib_mpp        ! MPP library
[1129]21   USE prtctl         ! Print control
[643]22
23   IMPLICIT NONE
24   PRIVATE
25
[1566]26   PUBLIC   dyn_adv_cen2   ! routine called by step.F90
[643]27
[3211]28   !! * Control permutation of array indices
29#  include "oce_ftrans.h90"
30#  include "dom_oce_ftrans.h90"
31
[643]32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
[2715]36   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]37   !! $Id$
[2715]38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[643]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
[1566]47      !!              and the general trend of the momentum equation.
[643]48      !!
49      !! ** Method  :   Trend evaluated using now fields (centered in time)
50      !!
[1566]51      !! ** Action  :   (ua,va) updated with the now vorticity term trend
[643]52      !!----------------------------------------------------------------------
[2715]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 
[3211]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
[2715]68      !
[1566]69      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
[2715]70      !
[1566]71      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[2715]72      REAL(wp) ::   zbu, zbv     ! local scalars
[643]73      !!----------------------------------------------------------------------
74
[2715]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,*) '~~~~~~~~~~~~'
[643]79      ENDIF
80
[2715]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
[1129]86      IF( l_trddyn ) THEN           ! Save ua and va trends
87         zfu_uw(:,:,:) = ua(:,:,:)
88         zfv_vw(:,:,:) = va(:,:,:)
89      ENDIF
[643]90
[1566]91      !                                      ! ====================== !
92      !                                      !  Horizontal advection  !
93      DO jk = 1, jpkm1                       ! ====================== !
94         !                                         ! horizontal volume fluxes
[643]95         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
96         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]97         !
98         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]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
[1566]106         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]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)
[1566]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
[643]115            END DO
116         END DO
[1566]117      END DO
118      !
119      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]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
[1566]126      !
[1129]127
[1566]128      !                                      ! ==================== !
129      !                                      !  Vertical advection  !
130      DO jk = 1, jpkm1                       ! ==================== !
131         !                                         ! Vertical volume fluxesÊ
[643]132         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]133         !
134         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
135            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]136            zfv_vw(:,:,jpk) = 0.e0
[1566]137            !                                           ! Surface value :
138            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]139               zfu_uw(:,:, 1 ) = 0.e0   
140               zfv_vw(:,:, 1 ) = 0.e0
[1566]141            ELSE                                             ! constant volume : advection through the surface
[643]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
[1566]149         ELSE                                      ! interior fluxes
[643]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
[1566]158      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]159         DO jj = 2, jpjm1 
160            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]161               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]162                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]163               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]164                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
165            END DO
166         END DO
167      END DO
[1566]168      !
169      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]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
[1566]174      !                                            ! Control print
[1129]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      !
[2715]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      !
[643]180   END SUBROUTINE dyn_adv_cen2
181
182   !!==============================================================================
183END MODULE dynadv_cen2
Note: See TracBrowser for help on using the repository browser.