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.
dynvor.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 41.5 KB
RevLine 
[3]1MODULE dynvor
2   !!======================================================================
3   !!                       ***  MODULE  dynvor  ***
4   !! Ocean dynamics: Update the momentum trend with the relative and
5   !!                 planetary vorticity trends
6   !!======================================================================
[2715]7   !! History :  OPA  ! 1989-12  (P. Andrich)  vor_ens: Original code
8   !!            5.0  ! 1991-11  (G. Madec) vor_ene, vor_mix: Original code
9   !!            6.0  ! 1996-01  (G. Madec)  s-coord, suppress work arrays
10   !!   NEMO     0.5  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-02  (G. Madec)  vor_een: Original code
12   !!             -   ! 2003-08  (G. Madec)  add vor_ctl
13   !!             -   ! 2005-11  (G. Madec)  add dyn_vor (new step architecture)
14   !!            2.0  ! 2006-11  (G. Madec)  flux form advection: add metric term
15   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
16   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[4990]17   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
[503]18   !!----------------------------------------------------------------------
[3]19
20   !!----------------------------------------------------------------------
[2528]21   !!   dyn_vor      : Update the momentum trend with the vorticity trend
22   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
23   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
24   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
25   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
26   !!   dyn_vor_init : set and control of the different vorticity option
[3]27   !!----------------------------------------------------------------------
[503]28   USE oce            ! ocean dynamics and tracers
29   USE dom_oce        ! ocean space and time domain
[3294]30   USE dommsk         ! ocean mask
[643]31   USE dynadv         ! momentum advection (use ln_dynadv_vec value)
[4990]32   USE trd_oce        ! trends: ocean variables
33   USE trddyn         ! trend manager: dynamics
[503]34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
35   USE prtctl         ! Print control
36   USE in_out_manager ! I/O manager
[3294]37   USE lib_mpp        ! MPP library
38   USE wrk_nemo       ! Memory Allocation
39   USE timing         ! Timing
[8280]40   USE lib_fortran
[3]41
[3294]42
[3]43   IMPLICIT NONE
44   PRIVATE
45
[2528]46   PUBLIC   dyn_vor        ! routine called by step.F90
47   PUBLIC   dyn_vor_init   ! routine called by opa.F90
[3]48
[4147]49   !                                   !!* Namelist namdyn_vor: vorticity term
50   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme
51   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme
52   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme
53   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme
[5029]54   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
[3]55
[503]56   INTEGER ::   nvor = 0   ! type of vorticity trend used
[643]57   INTEGER ::   ncor = 1   ! coriolis
58   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
59   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
[455]60
[3]61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
[2528]65   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]66   !! $Id$
[2715]67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]68   !!----------------------------------------------------------------------
69CONTAINS
70
[455]71   SUBROUTINE dyn_vor( kt )
[3]72      !!----------------------------------------------------------------------
73      !!
[455]74      !! ** Purpose :   compute the lateral ocean tracer physics.
75      !!
76      !! ** Action : - Update (ua,va) with the now vorticity term trend
[503]77      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]78      !!               and planetary vorticity trends) and send them to trd_dyn
79      !!               for futher diagnostics (l_trddyn=T)
[503]80      !!----------------------------------------------------------------------
[3294]81      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
[2715]82      !
[3294]83      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
[455]84      !!----------------------------------------------------------------------
[2715]85      !
[3294]86      IF( nn_timing == 1 )  CALL timing_start('dyn_vor')
87      !
88      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
89      !
[643]90      !                                          ! vorticity term
[455]91      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
[643]92      !
[455]93      CASE ( -1 )                                      ! esopa: test all possibility with control print
[643]94         CALL vor_ene( kt, ntot, ua, va )
[503]95         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
96            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]97         CALL vor_ens( kt, ntot, ua, va )
[503]98         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
99            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[455]100         CALL vor_mix( kt )
[503]101         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
102            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]103         CALL vor_een( kt, ntot, ua, va )
[503]104         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
105            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]106         !
[455]107      CASE ( 0 )                                       ! energy conserving scheme
108         IF( l_trddyn )   THEN
109            ztrdu(:,:,:) = ua(:,:,:)
110            ztrdv(:,:,:) = va(:,:,:)
[643]111            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
[455]112            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
113            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]114            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]115            ztrdu(:,:,:) = ua(:,:,:)
116            ztrdv(:,:,:) = va(:,:,:)
[643]117            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
[455]118            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
119            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]120            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]121         ELSE
[643]122            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
[455]123         ENDIF
[643]124         !
[455]125      CASE ( 1 )                                       ! enstrophy conserving scheme
126         IF( l_trddyn )   THEN   
127            ztrdu(:,:,:) = ua(:,:,:)
128            ztrdv(:,:,:) = va(:,:,:)
[643]129            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
[455]130            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
131            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]132            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]133            ztrdu(:,:,:) = ua(:,:,:)
134            ztrdv(:,:,:) = va(:,:,:)
[643]135            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
[455]136            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
137            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]138            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]139         ELSE
[643]140            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
[455]141         ENDIF
[643]142         !
[455]143      CASE ( 2 )                                       ! mixed ene-ens scheme
144         IF( l_trddyn )   THEN
145            ztrdu(:,:,:) = ua(:,:,:)
146            ztrdv(:,:,:) = va(:,:,:)
[643]147            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
[455]148            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
149            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]150            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]151            ztrdu(:,:,:) = ua(:,:,:)
152            ztrdv(:,:,:) = va(:,:,:)
[643]153            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
[455]154            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
155            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]156            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]157         ELSE
158            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
159         ENDIF
[643]160         !
[455]161      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
162         IF( l_trddyn )   THEN
163            ztrdu(:,:,:) = ua(:,:,:)
164            ztrdv(:,:,:) = va(:,:,:)
[643]165            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
[455]166            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
167            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]168            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]169            ztrdu(:,:,:) = ua(:,:,:)
170            ztrdv(:,:,:) = va(:,:,:)
[643]171            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
[455]172            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
173            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]174            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]175         ELSE
[643]176            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
[455]177         ENDIF
[643]178         !
[455]179      END SELECT
[2715]180      !
[455]181      !                       ! print sum trends (used for debugging)
[2715]182      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
[455]183         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]184      !
[3294]185      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
186      !
187      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor')
188      !
[455]189   END SUBROUTINE dyn_vor
190
191
[643]192   SUBROUTINE vor_ene( kt, kvor, pua, pva )
[455]193      !!----------------------------------------------------------------------
194      !!                  ***  ROUTINE vor_ene  ***
195      !!
[3]196      !! ** Purpose :   Compute the now total vorticity trend and add it to
197      !!      the general trend of the momentum equation.
198      !!
199      !! ** Method  :   Trend evaluated using now fields (centered in time)
200      !!      and the Sadourny (1975) flux form formulation : conserves the
201      !!      horizontal kinetic energy.
202      !!      The trend of the vorticity term is given by:
[455]203      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]204      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
205      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
206      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
207      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
208      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
209      !!      Add this trend to the general momentum trend (ua,va):
210      !!          (ua,va) = (ua,va) + ( voru , vorv )
211      !!
212      !! ** Action : - Update (ua,va) with the now vorticity term trend
213      !!
[503]214      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]215      !!----------------------------------------------------------------------
[2715]216      !
[643]217      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
218      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
[1438]219      !                                                           ! =nrvm (relative vorticity or metric)
[643]220      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
221      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
[2715]222      !
223      INTEGER  ::   ji, jj, jk   ! dummy loop indices
224      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
[3294]225      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
[3]226      !!----------------------------------------------------------------------
[3294]227      !
228      IF( nn_timing == 1 )  CALL timing_start('vor_ene')
229      !
230      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
231      !
[52]232      IF( kt == nit000 ) THEN
233         IF(lwp) WRITE(numout,*)
[455]234         IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
235         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[11101]236         IF(lflush) CALL flush(numout)
[52]237      ENDIF
[3]238
[1438]239      zfact2 = 0.5 * 0.5      ! Local constant initialization
[216]240
[455]241!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]242      !                                                ! ===============
243      DO jk = 1, jpkm1                                 ! Horizontal slab
244         !                                             ! ===============
[1438]245         !
[3]246         ! Potential vorticity and horizontal fluxes
247         ! -----------------------------------------
[643]248         SELECT CASE( kvor )      ! vorticity considered
249         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
250         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
251         CASE ( 3 )                                                ! metric term
252            DO jj = 1, jpjm1
253               DO ji = 1, fs_jpim1   ! vector opt.
254                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
255                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
256                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
257               END DO
258            END DO
259         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
260         CASE ( 5 )                                                ! total (coriolis + metric)
261            DO jj = 1, jpjm1
262               DO ji = 1, fs_jpim1   ! vector opt.
263                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
264                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
265                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
266                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
267                       &       )
268               END DO
269            END DO
[455]270         END SELECT
271
272         IF( ln_sco ) THEN
273            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
[3]274            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
275            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
276         ELSE
277            zwx(:,:) = e2u(:,:) * un(:,:,jk)
278            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
279         ENDIF
280
281         ! Compute and add the vorticity term trend
282         ! ----------------------------------------
283         DO jj = 2, jpjm1
284            DO ji = fs_2, fs_jpim1   ! vector opt.
285               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
286               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  )
287               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
288               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
[455]289               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
290               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
[3]291            END DO 
292         END DO 
293         !                                             ! ===============
294      END DO                                           !   End of slab
295      !                                                ! ===============
[3294]296      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
[2715]297      !
[3294]298      IF( nn_timing == 1 )  CALL timing_stop('vor_ene')
299      !
[455]300   END SUBROUTINE vor_ene
[216]301
302
[455]303   SUBROUTINE vor_mix( kt )
[3]304      !!----------------------------------------------------------------------
[455]305      !!                 ***  ROUTINE vor_mix  ***
[3]306      !!
307      !! ** Purpose :   Compute the now total vorticity trend and add it to
308      !!      the general trend of the momentum equation.
309      !!
310      !! ** Method  :   Trend evaluated using now fields (centered in time)
311      !!      Mixte formulation : conserves the potential enstrophy of a hori-
312      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
313      !!      ticity term and the horizontal kinetic energy for (f x uh), the
314      !!      coriolis term. the now trend of the vorticity term is given by:
[455]315      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]316      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
317      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
318      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
319      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
320      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
321      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
322      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
323      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
324      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
325      !!      Add this now trend to the general momentum trend (ua,va):
326      !!          (ua,va) = (ua,va) + ( voru , vorv )
327      !!
328      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
329      !!
[503]330      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]331      !!----------------------------------------------------------------------
[2715]332      !
[503]333      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
[2715]334      !
[1438]335      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[2715]336      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
337      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
[3294]338      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[3]339      !!----------------------------------------------------------------------
[3294]340      !
341      IF( nn_timing == 1 )  CALL timing_start('vor_mix')
342      !
343      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) 
344      !
[52]345      IF( kt == nit000 ) THEN
346         IF(lwp) WRITE(numout,*)
[455]347         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
348         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[11101]349         IF(lflush) CALL flush(numout)
[52]350      ENDIF
[3]351
[1438]352      zfact1 = 0.5 * 0.25      ! Local constant initialization
[3]353      zfact2 = 0.5 * 0.5
354
[455]355!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
[3]356      !                                                ! ===============
357      DO jk = 1, jpkm1                                 ! Horizontal slab
358         !                                             ! ===============
[1438]359         !
[3]360         ! Relative and planetary potential vorticity and horizontal fluxes
361         ! ----------------------------------------------------------------
[455]362         IF( ln_sco ) THEN       
[643]363            IF( ln_dynadv_vec ) THEN
364               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
365            ELSE                       
366               DO jj = 1, jpjm1
367                  DO ji = 1, fs_jpim1   ! vector opt.
368                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
369                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
370                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
371                  END DO
372               END DO
373            ENDIF
[3]374            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
375            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
376            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
377         ELSE
[643]378            IF( ln_dynadv_vec ) THEN
379               zww(:,:) = rotn(:,:,jk)
380            ELSE                       
381               DO jj = 1, jpjm1
382                  DO ji = 1, fs_jpim1   ! vector opt.
383                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
384                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
385                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
386                  END DO
387               END DO
388            ENDIF
389            zwz(:,:) = ff (:,:)
[3]390            zwx(:,:) = e2u(:,:) * un(:,:,jk)
391            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
392         ENDIF
393
394         ! Compute and add the vorticity term trend
395         ! ----------------------------------------
396         DO jj = 2, jpjm1
397            DO ji = fs_2, fs_jpim1   ! vector opt.
398               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
399               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
400               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
401               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
402               ! enstrophy conserving formulation for relative vorticity term
403               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
404               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
405               ! energy conserving formulation for planetary vorticity term
406               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
407               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
[503]408               ! mixed vorticity trend added to the momentum trends
[3]409               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
410               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
411            END DO 
412         END DO 
413         !                                             ! ===============
414      END DO                                           !   End of slab
415      !                                                ! ===============
[3294]416      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) 
[2715]417      !
[3294]418      IF( nn_timing == 1 )  CALL timing_stop('vor_mix')
419      !
[455]420   END SUBROUTINE vor_mix
[216]421
422
[643]423   SUBROUTINE vor_ens( kt, kvor, pua, pva )
[3]424      !!----------------------------------------------------------------------
[455]425      !!                ***  ROUTINE vor_ens  ***
[3]426      !!
427      !! ** Purpose :   Compute the now total vorticity trend and add it to
428      !!      the general trend of the momentum equation.
429      !!
430      !! ** Method  :   Trend evaluated using now fields (centered in time)
431      !!      and the Sadourny (1975) flux FORM formulation : conserves the
432      !!      potential enstrophy of a horizontally non-divergent flow. the
433      !!      trend of the vorticity term is given by:
[455]434      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
[3]435      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
436      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
437      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
438      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
439      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
440      !!      Add this trend to the general momentum trend (ua,va):
441      !!          (ua,va) = (ua,va) + ( voru , vorv )
442      !!
443      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
444      !!
[503]445      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]446      !!----------------------------------------------------------------------
[2715]447      !
[643]448      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
449      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
450         !                                                        ! =nrvm (relative vorticity or metric)
451      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
452      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
[2715]453      !
[503]454      INTEGER  ::   ji, jj, jk           ! dummy loop indices
455      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
[3294]456      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[3]457      !!----------------------------------------------------------------------
[3294]458      !
459      IF( nn_timing == 1 )  CALL timing_start('vor_ens')
460      !
461      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) 
462      !
[52]463      IF( kt == nit000 ) THEN
464         IF(lwp) WRITE(numout,*)
[455]465         IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
466         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[11101]467         IF(lflush) CALL flush(numout)
[52]468      ENDIF
[3]469
[1438]470      zfact1 = 0.5 * 0.25      ! Local constant initialization
[3]471
[455]472!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]473      !                                                ! ===============
474      DO jk = 1, jpkm1                                 ! Horizontal slab
475         !                                             ! ===============
[1438]476         !
[3]477         ! Potential vorticity and horizontal fluxes
478         ! -----------------------------------------
[643]479         SELECT CASE( kvor )      ! vorticity considered
480         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
481         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
482         CASE ( 3 )                                                ! metric term
483            DO jj = 1, jpjm1
484               DO ji = 1, fs_jpim1   ! vector opt.
485                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
486                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
487                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
488               END DO
489            END DO
490         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
491         CASE ( 5 )                                                ! total (coriolis + metric)
492            DO jj = 1, jpjm1
493               DO ji = 1, fs_jpim1   ! vector opt.
494                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
495                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
496                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
[1438]497                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
[643]498                       &       )
499               END DO
500            END DO
[455]501         END SELECT
[1438]502         !
[455]503         IF( ln_sco ) THEN
[3]504            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
505               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
[455]506                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
507                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
508                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
[3]509               END DO
510            END DO
511         ELSE
512            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
513               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
[455]514                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
515                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
[3]516               END DO
517            END DO
518         ENDIF
[1438]519         !
[3]520         ! Compute and add the vorticity term trend
521         ! ----------------------------------------
522         DO jj = 2, jpjm1
523            DO ji = fs_2, fs_jpim1   ! vector opt.
[455]524               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
[503]525                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
[455]526               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
[503]527                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
[455]528               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
529               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
[3]530            END DO 
531         END DO 
532         !                                             ! ===============
533      END DO                                           !   End of slab
534      !                                                ! ===============
[3294]535      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) 
[2715]536      !
[3294]537      IF( nn_timing == 1 )  CALL timing_stop('vor_ens')
538      !
[455]539   END SUBROUTINE vor_ens
[216]540
541
[643]542   SUBROUTINE vor_een( kt, kvor, pua, pva )
[108]543      !!----------------------------------------------------------------------
[455]544      !!                ***  ROUTINE vor_een  ***
[108]545      !!
546      !! ** Purpose :   Compute the now total vorticity trend and add it to
547      !!      the general trend of the momentum equation.
548      !!
549      !! ** Method  :   Trend evaluated using now fields (centered in time)
[1438]550      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]551      !!      both the horizontal kinetic energy and the potential enstrophy
[1438]552      !!      when horizontal divergence is zero (see the NEMO documentation)
553      !!      Add this trend to the general momentum trend (ua,va).
[108]554      !!
555      !! ** Action : - Update (ua,va) with the now vorticity term trend
556      !!
[503]557      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
558      !!----------------------------------------------------------------------
[2715]559      !
[643]560      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
561      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
[1438]562      !                                                           ! =nrvm (relative vorticity or metric)
[643]563      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
564      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
[218]565      !!
[3294]566      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
567      INTEGER  ::   ierr                                          ! local integer
568      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
[4292]569      REAL(wp) ::   zmsk, ze3                                     ! local scalars
[3294]570      !                                                           !  3D workspace
571      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
572      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
573#if defined key_vvl
574      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
[4292]575#else
[3294]576      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
[1438]577#endif
[108]578      !!----------------------------------------------------------------------
[3294]579      !
580      IF( nn_timing == 1 )  CALL timing_start('vor_een')
581      !
582      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
583      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
584#if defined key_vvl
585      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
586#endif
587      !
[108]588      IF( kt == nit000 ) THEN
589         IF(lwp) WRITE(numout,*)
[455]590         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
591         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[11101]592         IF(lflush) CALL flush(numout)
[3802]593#if ! defined key_vvl
594         IF( .NOT.ALLOCATED(ze3f) ) THEN
[2715]595            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
596            IF( lk_mpp    )   CALL mpp_sum ( ierr )
597            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
598         ENDIF
[4990]599         ze3f(:,:,:) = 0._wp
[3802]600#endif
[1438]601      ENDIF
[108]602
[4292]603      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
[5029]604
605         IF( ln_dynvor_een_old ) THEN ! original formulation
606            DO jk = 1, jpk
607               DO jj = 1, jpjm1
608                  DO ji = 1, jpim1
609                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
610                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
611                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
612                  END DO
[108]613               END DO
614            END DO
[5029]615         ELSE ! new formulation from NEMO 3.6
616            DO jk = 1, jpk
617               DO jj = 1, jpjm1
618                  DO ji = 1, jpim1
619                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
620                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
621                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
622                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
623                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
624                  END DO
625               END DO
626            END DO
627         ENDIF
628
[108]629         CALL lbc_lnk( ze3f, 'F', 1. )
630      ENDIF
631
[2715]632      zfac12 = 1._wp / 12._wp    ! Local constant initialization
[216]633
[108]634     
[455]635!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
[108]636      !                                                ! ===============
637      DO jk = 1, jpkm1                                 ! Horizontal slab
638         !                                             ! ===============
639         
640         ! Potential vorticity and horizontal fluxes
641         ! -----------------------------------------
[643]642         SELECT CASE( kvor )      ! vorticity considered
[1438]643         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
644            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
645         CASE ( 2 )                                                ! relative  vorticity
646            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
[643]647         CASE ( 3 )                                                ! metric term
648            DO jj = 1, jpjm1
649               DO ji = 1, fs_jpim1   ! vector opt.
650                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
651                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
652                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
653               END DO
654            END DO
[1516]655            CALL lbc_lnk( zwz, 'F', 1. )
656        CASE ( 4 )                                                ! total (relative + planetary vorticity)
[1438]657            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
[643]658         CASE ( 5 )                                                ! total (coriolis + metric)
659            DO jj = 1, jpjm1
660               DO ji = 1, fs_jpim1   ! vector opt.
661                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
662                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
663                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
[1438]664                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
[643]665                       &       ) * ze3f(ji,jj,jk)
666               END DO
667            END DO
[1516]668            CALL lbc_lnk( zwz, 'F', 1. )
[455]669         END SELECT
670
[108]671         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
672         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
673
674         ! Compute and add the vorticity term trend
675         ! ----------------------------------------
[1438]676         jj = 2
677         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
[108]678         DO ji = 2, jpi   
679               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
680               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
681               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
682               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
683         END DO
684         DO jj = 3, jpj
[1694]685            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
[108]686               ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
687               ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
688               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
689               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
690            END DO
691         END DO
692         DO jj = 2, jpjm1
693            DO ji = fs_2, fs_jpim1   ! vector opt.
694               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
695                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
696               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
697                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
[455]698               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
699               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]700            END DO 
701         END DO 
702         !                                             ! ===============
703      END DO                                           !   End of slab
704      !                                                ! ===============
[3294]705      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
706      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
707#if defined key_vvl
708      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
709#endif
[2715]710      !
[3294]711      IF( nn_timing == 1 )  CALL timing_stop('vor_een')
712      !
[455]713   END SUBROUTINE vor_een
[216]714
715
[2528]716   SUBROUTINE dyn_vor_init
[3]717      !!---------------------------------------------------------------------
[2528]718      !!                  ***  ROUTINE dyn_vor_init  ***
[3]719      !!
720      !! ** Purpose :   Control the consistency between cpp options for
[1438]721      !!              tracer advection schemes
[3]722      !!----------------------------------------------------------------------
[2715]723      INTEGER ::   ioptio          ! local integer
[3294]724      INTEGER ::   ji, jj, jk      ! dummy loop indices
[4147]725      INTEGER ::   ios             ! Local integer output status for namelist read
[2715]726      !!
[5029]727      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
[3]728      !!----------------------------------------------------------------------
729
[4147]730      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
731      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
732901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
[3]733
[4147]734      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
735      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
736902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
[11101]737      IF(lwm .AND. nprint >2) WRITE ( numond, namdyn_vor )
[4147]738
[503]739      IF(lwp) THEN                    ! Namelist print
[3]740         WRITE(numout,*)
[2528]741         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
742         WRITE(numout,*) '~~~~~~~~~~~~'
[4147]743         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
[503]744         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
745         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
746         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
747         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
[5029]748         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
[11101]749         IF(lflush) CALL flush(numout)
[52]750      ENDIF
751
[3294]752      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
753      ! at angles with three ocean points and one land point
754      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
755         DO jk = 1, jpk
756            DO jj = 2, jpjm1
757               DO ji = 2, jpim1
758                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
759                      fmask(ji,jj,jk) = 1._wp
760               END DO
761            END DO
762         END DO
763          !
764          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
765          !
766      ENDIF
767
[503]768      ioptio = 0                     ! Control of vorticity scheme options
769      IF( ln_dynvor_ene )   ioptio = ioptio + 1
770      IF( ln_dynvor_ens )   ioptio = ioptio + 1
771      IF( ln_dynvor_mix )   ioptio = ioptio + 1
772      IF( ln_dynvor_een )   ioptio = ioptio + 1
[5029]773      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
[503]774      IF( lk_esopa      )   ioptio =          1
775
776      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
777
[643]778      !                              ! Set nvor (type of scheme for vorticity)
[503]779      IF( ln_dynvor_ene )   nvor =  0
780      IF( ln_dynvor_ens )   nvor =  1
781      IF( ln_dynvor_mix )   nvor =  2
[5029]782      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
[503]783      IF( lk_esopa      )   nvor = -1
784     
[643]785      !                              ! Set ncor, nrvm, ntot (type of vorticity)
786      IF(lwp) WRITE(numout,*)
787      ncor = 1
788      IF( ln_dynadv_vec ) THEN     
789         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
790         nrvm = 2
791         ntot = 4
792      ELSE                       
793         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
794         nrvm = 3
795         ntot = 5
796      ENDIF
797     
[503]798      IF(lwp) THEN                   ! Print the choice
799         WRITE(numout,*)
[643]800         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
801         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
802         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
803         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
[503]804         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
[11101]805         IF(lflush) CALL flush(numout)
[3]806      ENDIF
[503]807      !
[2528]808   END SUBROUTINE dyn_vor_init
[3]809
[503]810   !!==============================================================================
[3]811END MODULE dynvor
Note: See TracBrowser for help on using the repository browser.