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.
dynhpg.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/dynhpg.F90 @ 3211

Last change on this file since 3211 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: 56.3 KB
Line 
1MODULE dynhpg
2   !!======================================================================
3   !!                       ***  MODULE  dynhpg  ***
4   !! Ocean dynamics:  hydrostatic pressure gradient trend
5   !!======================================================================
6   !! History :  OPA  !  1987-09  (P. Andrich, M.-A. Foujols)  hpg_zco: Original code
7   !!            5.0  !  1991-11  (G. Madec)
8   !!            7.0  !  1996-01  (G. Madec)  hpg_sco: Original code for s-coordinates
9   !!            8.0  !  1997-05  (G. Madec)  split dynber into dynkeg and dynhpg
10   !!            8.5  !  2002-07  (G. Madec)  F90: Free form and module
11   !!            8.5  !  2002-08  (A. Bozec)  hpg_zps: Original code
12   !!   NEMO     1.0  !  2005-10  (A. Beckmann, B.W. An)  various s-coordinate options
13   !!                 !         Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot
14   !!             -   !  2005-11  (G. Madec) style & small optimisation
15   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   dyn_hpg      : update the momentum trend with the now horizontal
20   !!                  gradient of the hydrostatic pressure
21   !!   dyn_hpg_init : initialisation and control of options
22   !!       hpg_zco  : z-coordinate scheme
23   !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
24   !!       hpg_sco  : s-coordinate (standard jacobian formulation)
25   !!       hpg_hel  : s-coordinate (helsinki modification)
26   !!       hpg_wdj  : s-coordinate (weighted density jacobian)
27   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
28   !!       hpg_rot  : s-coordinate (ROTated axes scheme)
29   !!----------------------------------------------------------------------
30   USE oce             ! ocean dynamics and tracers
31   USE dom_oce         ! ocean space and time domain
32   USE phycst          ! physical constants
33   USE trdmod          ! ocean dynamics trends
34   USE trdmod_oce      ! ocean variables trends
35   USE in_out_manager  ! I/O manager
36   USE prtctl          ! Print control
37   USE lbclnk          ! lateral boundary condition
38   USE lib_mpp         ! MPP library
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   dyn_hpg        ! routine called by step module
44   PUBLIC   dyn_hpg_init   ! routine called by opa module
45
46   !                                              !!* Namelist namdyn_hpg : hydrostatic pressure gradient
47   LOGICAL , PUBLIC ::   ln_hpg_zco    = .TRUE.    !: z-coordinate - full steps
48   LOGICAL , PUBLIC ::   ln_hpg_zps    = .FALSE.   !: z-coordinate - partial steps (interpolation)
49   LOGICAL , PUBLIC ::   ln_hpg_sco    = .FALSE.   !: s-coordinate (standard jacobian formulation)
50   LOGICAL , PUBLIC ::   ln_hpg_hel    = .FALSE.   !: s-coordinate (helsinki modification)
51   LOGICAL , PUBLIC ::   ln_hpg_wdj    = .FALSE.   !: s-coordinate (weighted density jacobian)
52   LOGICAL , PUBLIC ::   ln_hpg_djc    = .FALSE.   !: s-coordinate (Density Jacobian with Cubic polynomial)
53   LOGICAL , PUBLIC ::   ln_hpg_rot    = .FALSE.   !: s-coordinate (ROTated axes scheme)
54   REAL(wp), PUBLIC ::   rn_gamma      = 0._wp     !: weighting coefficient
55   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag
56
57   INTEGER  ::   nhpg  =  0   ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)
58
59   !! * Control permutation of array indices
60#  include "oce_ftrans.h90"
61#  include "dom_oce_ftrans.h90"
62
63   !! * Substitutions
64#  include "domzgr_substitute.h90"
65#  include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
68   !! $Id$
69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE dyn_hpg( kt )
74      !!---------------------------------------------------------------------
75      !!                  ***  ROUTINE dyn_hpg  ***
76      !!
77      !! ** Method  :   Call the hydrostatic pressure gradient routine
78      !!              using the scheme defined in the namelist
79      !!   
80      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
81      !!             - Save the trend (l_trddyn=T)
82      !!----------------------------------------------------------------------
83      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
84      USE wrk_nemo, ONLY:   ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2   ! 3D workspace
85      !! DCSE_NEMO: need additional directives for renamed module variables
86!FTRANS ztrdu :I :I :z
87!FTRANS ztrdv :I :I :z
88
89      !!
90      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
91      !!----------------------------------------------------------------------
92      !
93      IF( wrk_in_use(3, 1,2) ) THEN
94         CALL ctl_stop('dyn_hpg: requested workspace arrays are unavailable')   ;   RETURN
95      ENDIF
96      !
97      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
98         ztrdu(:,:,:) = ua(:,:,:) 
99         ztrdv(:,:,:) = va(:,:,:) 
100      ENDIF     
101      !
102      SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation
103      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
104      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
105      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
106      CASE (  3 )   ;   CALL hpg_hel    ( kt )      ! s-coordinate (helsinki modification)
107      CASE (  4 )   ;   CALL hpg_wdj    ( kt )      ! s-coordinate (weighted density jacobian)
108      CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
109      CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme)
110      END SELECT
111      !
112      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
113         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
114         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
115         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
116      ENDIF         
117      !
118      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
119         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
120      !
121      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_hpg: failed to release workspace arrays')
122      !
123   END SUBROUTINE dyn_hpg
124
125
126   SUBROUTINE dyn_hpg_init
127      !!----------------------------------------------------------------------
128      !!                 ***  ROUTINE dyn_hpg_init  ***
129      !!
130      !! ** Purpose :   initializations for the hydrostatic pressure gradient
131      !!              computation and consistency control
132      !!
133      !! ** Action  :   Read the namelist namdyn_hpg and check the consistency
134      !!      with the type of vertical coordinate used (zco, zps, sco)
135      !!----------------------------------------------------------------------
136      INTEGER ::   ioptio = 0      ! temporary integer
137      !!
138      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,    &
139         &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma  , ln_dynhpg_imp
140      !!----------------------------------------------------------------------
141      !
142      REWIND( numnam )               ! Read Namelist namdyn_hpg
143      READ  ( numnam, namdyn_hpg )
144      !
145      IF(lwp) THEN                   ! Control print
146         WRITE(numout,*)
147         WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
148         WRITE(numout,*) '~~~~~~~~~~~~'
149         WRITE(numout,*) '   Namelist namdyn_hpg : choice of hpg scheme'
150         WRITE(numout,*) '      z-coord. - full steps                             ln_hpg_zco    = ', ln_hpg_zco
151         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps
152         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco
153         WRITE(numout,*) '      s-coord. (helsinki modification)                  ln_hpg_hel    = ', ln_hpg_hel
154         WRITE(numout,*) '      s-coord. (weighted density jacobian)              ln_hpg_wdj    = ', ln_hpg_wdj
155         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
156         WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot
157         WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma
158         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp
159      ENDIF
160      !
161      IF( lk_vvl .AND. .NOT. ln_hpg_sco )   &
162         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco')
163      !
164      !                               ! Set nhpg from ln_hpg_... flags
165      IF( ln_hpg_zco )   nhpg = 0
166      IF( ln_hpg_zps )   nhpg = 1
167      IF( ln_hpg_sco )   nhpg = 2
168      IF( ln_hpg_hel )   nhpg = 3
169      IF( ln_hpg_wdj )   nhpg = 4
170      IF( ln_hpg_djc )   nhpg = 5
171      IF( ln_hpg_rot )   nhpg = 6
172      !
173      !                               ! Consitency check
174      ioptio = 0 
175      IF( ln_hpg_zco )   ioptio = ioptio + 1
176      IF( ln_hpg_zps )   ioptio = ioptio + 1
177      IF( ln_hpg_sco )   ioptio = ioptio + 1
178      IF( ln_hpg_hel )   ioptio = ioptio + 1
179      IF( ln_hpg_wdj )   ioptio = ioptio + 1
180      IF( ln_hpg_djc )   ioptio = ioptio + 1
181      IF( ln_hpg_rot )   ioptio = ioptio + 1
182      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
183      !
184   END SUBROUTINE dyn_hpg_init
185
186
187   SUBROUTINE hpg_zco( kt )
188      !!---------------------------------------------------------------------
189      !!                  ***  ROUTINE hpg_zco  ***
190      !!
191      !! ** Method  :   z-coordinate case, levels are horizontal surfaces.
192      !!      The now hydrostatic pressure gradient at a given level, jk,
193      !!      is computed by taking the vertical integral of the in-situ
194      !!      density gradient along the model level from the suface to that
195      !!      level:    zhpi = grav .....
196      !!                zhpj = grav .....
197      !!      add it to the general momentum trend (ua,va).
198      !!            ua = ua - 1/e1u * zhpi
199      !!            va = va - 1/e2v * zhpj
200      !!
201      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
202      !!----------------------------------------------------------------------
203      USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace
204      !! DCSE_NEMO: need additional directives for renamed module variables
205!FTRANS zhpi :I :I :z
206!FTRANS zhpj :I :I :z
207
208      !!
209      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
210      !!
211      INTEGER  ::   ji, jj, jk       ! dummy loop indices
212      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
213      !!----------------------------------------------------------------------
214     
215      IF( kt == nit000 ) THEN
216         IF(lwp) WRITE(numout,*)
217         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
218         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
219      ENDIF
220     
221      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
222
223      ! Surface value
224      DO jj = 2, jpjm1
225         DO ji = fs_2, fs_jpim1   ! vector opt.
226            zcoef1 = zcoef0 * fse3w(ji,jj,1)
227            ! hydrostatic pressure gradient
228            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
229            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
230            ! add to the general momentum trend
231            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
232            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
233         END DO
234      END DO
235      !
236      ! interior value (2=<jk=<jpkm1)
237#if defined key_z_first
238      DO jj = 2, jpjm1
239         DO ji = 2, jpim1
240            DO jk = 2, jpkm1
241#else
242      DO jk = 2, jpkm1
243         DO jj = 2, jpjm1
244            DO ji = fs_2, fs_jpim1   ! vector opt.
245#endif
246               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
247               ! hydrostatic pressure gradient
248               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
249                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
250                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
251
252               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
253                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
254                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
255               ! add to the general momentum trend
256               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
257               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
258            END DO
259         END DO
260      END DO
261      !
262   END SUBROUTINE hpg_zco
263
264
265   SUBROUTINE hpg_zps( kt )
266      !!---------------------------------------------------------------------
267      !!                 ***  ROUTINE hpg_zps  ***
268      !!                   
269      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
270      !!
271      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
272      !!----------------------------------------------------------------------
273      USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace
274      !! DCSE_NEMO: need additional directives for renamed module variables
275!FTRANS zhpi :I :I :z
276!FTRANS zhpj :I :I :z
277      !!
278      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
279      !!
280      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
281      INTEGER  ::   iku, ikv                         ! temporary integers
282      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
283      !!----------------------------------------------------------------------
284
285      IF( kt == nit000 ) THEN
286         IF(lwp) WRITE(numout,*)
287         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
288         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
289      ENDIF
290
291      ! Local constant initialization
292      zcoef0 = - grav * 0.5_wp
293
294      !  Surface value (also valid in partial step case)
295      DO jj = 2, jpjm1
296         DO ji = fs_2, fs_jpim1   ! vector opt.
297            zcoef1 = zcoef0 * fse3w(ji,jj,1)
298            ! hydrostatic pressure gradient
299            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
300            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
301            ! add to the general momentum trend
302            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
303            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
304         END DO
305      END DO
306
307      ! interior value (2=<jk=<jpkm1)
308#if defined key_z_first
309      DO jj = 2, jpjm1
310         DO ji = 2, jpim1
311            DO jk = 2, jpkm1
312#else
313      DO jk = 2, jpkm1
314         DO jj = 2, jpjm1
315            DO ji = fs_2, fs_jpim1   ! vector opt.
316#endif
317               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
318               ! hydrostatic pressure gradient
319               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
320                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
321                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
322
323               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
324                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
325                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
326               ! add to the general momentum trend
327               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
328               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
329            END DO
330         END DO
331      END DO
332
333      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
334# if defined key_vectopt_loop
335         jj = 1
336         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
337# else
338      DO jj = 2, jpjm1
339         DO ji = 2, jpim1
340# endif
341            iku = mbku(ji,jj)
342            ikv = mbkv(ji,jj)
343            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) )
344            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) )
345            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
346               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
347               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
348                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
349               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
350            ENDIF
351            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
352               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
353               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
354                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
355               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
356            ENDIF
357# if ! defined key_vectopt_loop
358         END DO
359# endif
360      END DO
361      !
362   END SUBROUTINE hpg_zps
363
364
365   SUBROUTINE hpg_sco( kt )
366      !!---------------------------------------------------------------------
367      !!                  ***  ROUTINE hpg_sco  ***
368      !!
369      !! ** Method  :   s-coordinate case. Jacobian scheme.
370      !!      The now hydrostatic pressure gradient at a given level, jk,
371      !!      is computed by taking the vertical integral of the in-situ
372      !!      density gradient along the model level from the suface to that
373      !!      level. s-coordinates (ln_sco): a corrective term is added
374      !!      to the horizontal pressure gradient :
375      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
376      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
377      !!      add it to the general momentum trend (ua,va).
378      !!         ua = ua - 1/e1u * zhpi
379      !!         va = va - 1/e2v * zhpj
380      !!
381      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
382      !!----------------------------------------------------------------------
383      USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace
384      !! DCSE_NEMO: need additional directives for renamed module variables
385!FTRANS zhpi :I :I :z
386!FTRANS zhpj :I :I :z
387      !!
388      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
389      !!
390      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
391      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
392      !!----------------------------------------------------------------------
393
394      IF( kt == nit000 ) THEN
395         IF(lwp) WRITE(numout,*)
396         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
397         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
398      ENDIF
399
400      ! Local constant initialization
401      zcoef0 = - grav * 0.5_wp
402      ! To use density and not density anomaly
403      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
404      ELSE                 ;     znad = 0._wp         ! Fixed volume
405      ENDIF
406
407      ! Surface value
408      DO jj = 2, jpjm1
409         DO ji = fs_2, fs_jpim1   ! vector opt.   
410            ! hydrostatic pressure gradient along s-surfaces
411            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
412               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
413            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
414               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
415            ! s-coordinate pressure gradient correction
416            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
417               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
418            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
419               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
420            ! add to the general momentum trend
421            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
422            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
423         END DO 
424      END DO   
425           
426      ! interior value (2=<jk=<jpkm1)
427#if defined key_z_first
428      DO jj = 2, jpjm1     
429         DO ji = 2, jpim1
430            DO jk = 2, jpkm1                                 
431#else
432      DO jk = 2, jpkm1                                 
433         DO jj = 2, jpjm1     
434            DO ji = fs_2, fs_jpim1   ! vector opt.     
435#endif
436               ! hydrostatic pressure gradient along s-surfaces
437               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
438                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
439                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
440               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
441                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
442                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
443               ! s-coordinate pressure gradient correction
444               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
445                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
446               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
447                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
448               ! add to the general momentum trend
449               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
450               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
451            END DO
452         END DO
453      END DO
454      !
455   END SUBROUTINE hpg_sco
456
457
458   SUBROUTINE hpg_hel( kt )
459      !!---------------------------------------------------------------------
460      !!                  ***  ROUTINE hpg_hel  ***
461      !!
462      !! ** Method  :   s-coordinate case.
463      !!      The now hydrostatic pressure gradient at a given level
464      !!      jk is computed by taking the vertical integral of the in-situ
465      !!      density gradient along the model level from the suface to that
466      !!      level. s-coordinates (ln_sco): a corrective term is added
467      !!      to the horizontal pressure gradient :
468      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
469      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
470      !!      add it to the general momentum trend (ua,va).
471      !!         ua = ua - 1/e1u * zhpi
472      !!         va = va - 1/e2v * zhpj
473      !!
474      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
475      !!             - Save the trend (l_trddyn=T)
476      !!----------------------------------------------------------------------
477      USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace
478      !! DCSE_NEMO: need additional directives for renamed module variables
479!FTRANS zhpi :I :I :z
480!FTRANS zhpj :I :I :z
481      !!
482      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
483      !!
484      INTEGER  ::   ji, jj, jk           ! dummy loop indices
485      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
486      !!----------------------------------------------------------------------
487
488      IF( kt == nit000 ) THEN
489         IF(lwp) WRITE(numout,*)
490         IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend'
491         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, helsinki modified scheme'
492      ENDIF
493
494      ! Local constant initialization
495      zcoef0 = - grav * 0.5_wp
496 
497      ! Surface value
498      DO jj = 2, jpjm1
499         DO ji = fs_2, fs_jpim1   ! vector opt.
500            ! hydrostatic pressure gradient along s-surfaces
501            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  &
502               &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) )
503            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)  &
504               &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) )
505            ! s-coordinate pressure gradient correction
506            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   &
507               &           * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
508            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   &
509               &           * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
510            ! add to the general momentum trend
511            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
512            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
513         END DO
514      END DO
515      !
516      ! interior value (2=<jk=<jpkm1)
517#if defined key_z_first
518      DO jj = 2, jpjm1
519         DO ji = 2, jpim1
520            DO jk = 2, jpkm1
521#else
522      DO jk = 2, jpkm1
523         DO jj = 2, jpjm1
524            DO ji = fs_2, fs_jpim1   ! vector opt.
525#endif
526               ! hydrostatic pressure gradient along s-surfaces
527               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
528                  &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk)     &
529                  &                                     -fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk)   ) &
530                  &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   &
531                  &                                     -fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1) )
532               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
533                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)     &
534                  &                                     -fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk)   ) &
535                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   &
536                  &                                     -fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1) )
537               ! s-coordinate pressure gradient correction
538               zuap = - zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )   &
539                  &            * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
540               zvap = - zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )   &
541                  &            * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
542               ! add to the general momentum trend
543               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
544               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
545            END DO
546         END DO
547      END DO
548      !
549   END SUBROUTINE hpg_hel
550
551
552   SUBROUTINE hpg_wdj( kt )
553      !!---------------------------------------------------------------------
554      !!                  ***  ROUTINE hpg_wdj  ***
555      !!
556      !! ** Method  :   Weighted Density Jacobian (wdj) scheme (song 1998)
557      !!      The weighting coefficients from the namelist parameter rn_gamma
558      !!      (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma
559      !!
560      !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.
561      !!----------------------------------------------------------------------
562      USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace
563      !! DCSE_NEMO: need additional directives for renamed module variables
564!FTRANS zhpi :I :I :z
565!FTRANS zhpj :I :I :z
566      !!
567      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
568      !!
569      INTEGER  ::   ji, jj, jk           ! dummy loop indices
570      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
571      REAL(wp) ::   zalph , zbeta        !    "         "
572      !!----------------------------------------------------------------------
573
574      IF( kt == nit000 ) THEN
575         IF(lwp) WRITE(numout,*)
576         IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend'
577         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   Weighted Density Jacobian'
578      ENDIF
579
580      ! Local constant initialization
581      zcoef0 = - grav * 0.5_wp
582      zalph  = 0.5_wp - rn_gamma    ! weighting coefficients (alpha=0.5-rn_gamma
583      zbeta  = 0.5_wp + rn_gamma    !                        (beta =1-alpha=0.5+rn_gamma
584
585      ! Surface value (no ponderation)
586      DO jj = 2, jpjm1
587         DO ji = fs_2, fs_jpim1   ! vector opt.
588            ! hydrostatic pressure gradient along s-surfaces
589            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3w(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)   &
590               &                                   - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  )
591            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3w(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   &
592               &                                   - fse3w(ji  ,jj  ,1) * rhd(ji,  jj  ,1)  )
593            ! s-coordinate pressure gradient correction
594            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   &
595               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
596            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   &
597               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
598            ! add to the general momentum trend
599            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
600            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
601         END DO
602      END DO
603
604      ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)
605#if defined key_z_first
606      DO jj = 2, jpjm1
607         DO ji = 2, jpim1
608            DO jk = 2, jpkm1
609#else
610      DO jk = 2, jpkm1
611         DO jj = 2, jpjm1
612            DO ji = fs_2, fs_jpim1   ! vector opt.
613#endif
614               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)                            &
615                  &           * (   (            fsde3w(ji+1,jj,jk  ) + fsde3w(ji,jj,jk  )        &
616                  &                            - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1)    )   &
617                  &               * (  zalph * ( rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1) )      &
618                  &                  + zbeta * ( rhd   (ji+1,jj,jk  ) - rhd   (ji,jj,jk  ) )  )   &
619                  &             -   (            rhd   (ji+1,jj,jk  ) + rhd   (ji,jj,jk  )        &
620                  &                           - rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1)     )   &
621                  &               * (  zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) )      &
622                  &                  + zbeta * ( fsde3w(ji+1,jj,jk  ) - fsde3w(ji,jj,jk  ) )  )  )
623               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)                            &
624                  &           * (   (           fsde3w(ji,jj+1,jk  ) + fsde3w(ji,jj,jk  )         &
625                  &                           - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1)     )   &
626                  &               * (  zalph * ( rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1) )      &
627                  &                  + zbeta * ( rhd   (ji,jj+1,jk  ) - rhd   (ji,jj,jk  ) )  )   &
628                  &             -   (            rhd   (ji,jj+1,jk  ) + rhd   (ji,jj,jk  )        &
629                  &                            - rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1)    )   &
630                  &               * (  zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) )      &
631                  &                  + zbeta * ( fsde3w(ji,jj+1,jk  ) - fsde3w(ji,jj,jk  ) )  )  )
632               ! add to the general momentum trend
633               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
634               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
635            END DO
636         END DO
637      END DO
638      !
639   END SUBROUTINE hpg_wdj
640
641
642   SUBROUTINE hpg_djc( kt )
643      !!---------------------------------------------------------------------
644      !!                  ***  ROUTINE hpg_djc  ***
645      !!
646      !! ** Method  :   Density Jacobian with Cubic polynomial scheme
647      !!
648      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
649      !!----------------------------------------------------------------------
650      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
651      USE oce     , ONLY:   zhpi  => ta        , zhpj => sa       ! (ta,sa) used as 3D workspace
652      USE wrk_nemo, ONLY:   drhox => wrk_3d_1  , dzx  => wrk_3d_2
653      USE wrk_nemo, ONLY:   drhou => wrk_3d_3  , dzu  => wrk_3d_4 , rho_i => wrk_3d_5
654      USE wrk_nemo, ONLY:   drhoy => wrk_3d_6  , dzy  => wrk_3d_7
655      USE wrk_nemo, ONLY:   drhov => wrk_3d_8  , dzv  => wrk_3d_9 , rho_j => wrk_3d_10
656      USE wrk_nemo, ONLY:   drhoz => wrk_3d_11 , dzz  => wrk_3d_12 
657      USE wrk_nemo, ONLY:   drhow => wrk_3d_13 , dzw  => wrk_3d_14
658      USE wrk_nemo, ONLY:   rho_k => wrk_3d_15
659      !! DCSE_NEMO: need additional directives for renamed module variables
660!FTRANS zhpi :I :I :z
661!FTRANS zhpj :I :I :z
662!FTRANS drhox :I :I :z
663!FTRANS dzx :I :I :z
664!FTRANS drhou :I :I :z
665!FTRANS dzu :I :I :z
666!FTRANS rho_i :I :I :z
667!FTRANS drhoy :I :I :z
668!FTRANS dzy :I :I :z
669!FTRANS drhov :I :I :z
670!FTRANS dzv :I :I :z
671!FTRANS rho_j :I :I :z
672!FTRANS drhoz :I :I :z
673!FTRANS dzz :I :I :z
674!FTRANS drhow :I :I :z
675!FTRANS dzw :I :I :z
676!FTRANS rho_k :I :I :z
677      !!
678      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
679      !!
680      INTEGER  ::   ji, jj, jk          ! dummy loop indices
681      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
682      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
683      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
684      !!----------------------------------------------------------------------
685
686      IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) THEN
687         CALL ctl_stop('dyn:hpg_djc: requested workspace arrays unavailable')   ;   RETURN
688      ENDIF
689
690      IF( kt == nit000 ) THEN
691         IF(lwp) WRITE(numout,*)
692         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
693         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
694      ENDIF
695
696      ! Local constant initialization
697      zcoef0 = - grav * 0.5_wp
698      z1_10  = 1._wp / 10._wp
699      z1_12  = 1._wp / 12._wp
700
701      !----------------------------------------------------------------------------------------
702      !  compute and store in provisional arrays elementary vertical and horizontal differences
703      !----------------------------------------------------------------------------------------
704
705!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
706
707#if defined key_z_first
708      DO jj = 2, jpjm1
709         DO ji = 2, jpim1
710            DO jk = 2, jpkm1
711#else
712      DO jk = 2, jpkm1
713         DO jj = 2, jpjm1
714            DO ji = fs_2, fs_jpim1   ! vector opt.
715#endif
716               drhoz(ji,jj,jk) = rhd   (ji  ,jj  ,jk) - rhd   (ji,jj,jk-1)
717               dzz  (ji,jj,jk) = fsde3w(ji  ,jj  ,jk) - fsde3w(ji,jj,jk-1)
718               drhox(ji,jj,jk) = rhd   (ji+1,jj  ,jk) - rhd   (ji,jj,jk  )
719               dzx  (ji,jj,jk) = fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk  )
720               drhoy(ji,jj,jk) = rhd   (ji  ,jj+1,jk) - rhd   (ji,jj,jk  )
721               dzy  (ji,jj,jk) = fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk  )
722            END DO
723         END DO
724      END DO
725
726      !-------------------------------------------------------------------------
727      ! compute harmonic averages using eq. 5.18
728      !-------------------------------------------------------------------------
729      zep = 1.e-15
730
731!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
732!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
733
734#if defined key_z_first
735      DO jj = 2, jpjm1
736         DO ji = 2, jpim1
737            DO jk = 2, jpkm1
738#else
739      DO jk = 2, jpkm1
740         DO jj = 2, jpjm1
741            DO ji = fs_2, fs_jpim1   ! vector opt.
742#endif
743               cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
744
745               cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  )
746               cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  )
747 
748               cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  )
749               cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  )
750
751               IF( cffw > zep) THEN
752                  drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   &
753                     &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
754               ELSE
755                  drhow(ji,jj,jk) = 0._wp
756               ENDIF
757
758               dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   &
759                  &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
760
761               IF( cffu > zep ) THEN
762                  drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   &
763                     &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
764               ELSE
765                  drhou(ji,jj,jk ) = 0._wp
766               ENDIF
767
768               IF( cffx > zep ) THEN
769                  dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   &
770                     &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
771               ELSE
772                  dzu(ji,jj,jk) = 0._wp
773               ENDIF
774
775               IF( cffv > zep ) THEN
776                  drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   &
777                     &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
778               ELSE
779                  drhov(ji,jj,jk) = 0._wp
780               ENDIF
781
782               IF( cffy > zep ) THEN
783                  dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   &
784                     &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
785               ELSE
786                  dzv(ji,jj,jk) = 0._wp
787               ENDIF
788
789            END DO
790         END DO
791      END DO
792
793      !----------------------------------------------------------------------------------
794      ! apply boundary conditions at top and bottom using 5.36-5.37
795      !----------------------------------------------------------------------------------
796      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  )
797      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  )
798      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  )
799
800      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
801      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
802      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
803
804
805      !--------------------------------------------------------------
806      ! Upper half of top-most grid box, compute and store
807      !-------------------------------------------------------------
808
809!!bug gm   :  e3w-de3w = 0.5*e3w  ....  and de3w(2)-de3w(1)=e3w(2) ....   to be verified
810!          true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
811
812      DO jj = 2, jpjm1
813         DO ji = fs_2, fs_jpim1   ! vector opt.
814            rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )               &
815               &                   * (  rhd(ji,jj,1)                                    &
816               &                     + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         &
817               &                              * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   &
818               &                              / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  ) 
819         END DO
820      END DO
821
822!!bug gm    : here also, simplification is possible
823!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
824
825#if defined key_z_first
826      DO jj = 2, jpjm1
827         DO ji = 2, jpim1
828            DO jk = 2, jpkm1
829#else
830      DO jk = 2, jpkm1
831         DO jj = 2, jpjm1
832            DO ji = fs_2, fs_jpim1   ! vector opt.
833#endif
834
835               rho_k(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj,jk) + rhd   (ji,jj,jk-1) )                                   &
836                  &                     * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) )                                   &
837                  &            - grav * z1_10 * (                                                                     &
838                  &     ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) )                                                     &
839                  &   * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
840                  &   - ( dzw   (ji,jj,jk) - dzw   (ji,jj,jk-1) )                                                     &
841                  &   * ( rhd   (ji,jj,jk) - rhd   (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   &
842                  &                             )
843
844               rho_i(ji,jj,jk) = zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )                                   &
845                  &                     * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) )                                   &
846                  &            - grav* z1_10 * (                                                                      &
847                  &     ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) )                                                     &
848                  &   * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
849                  &   - ( dzu   (ji+1,jj,jk) - dzu   (ji,jj,jk) )                                                     &
850                  &   * ( rhd   (ji+1,jj,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   &
851                  &                            )
852
853               rho_j(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )                                   &
854                  &                     * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) )                                   &
855                  &            - grav* z1_10 * (                                                                      &
856                  &     ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) )                                                     &
857                  &   * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
858                  &   - ( dzv   (ji,jj+1,jk) - dzv   (ji,jj,jk) )                                                     &
859                  &   * ( rhd   (ji,jj+1,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   &
860                  &                            )
861
862            END DO
863         END DO
864      END DO
865      CALL lbc_lnk(rho_k,'W',1.)
866      CALL lbc_lnk(rho_i,'U',1.)
867      CALL lbc_lnk(rho_j,'V',1.)
868
869
870      ! ---------------
871      !  Surface value
872      ! ---------------
873      DO jj = 2, jpjm1
874         DO ji = fs_2, fs_jpim1   ! vector opt.
875            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
876            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
877            ! add to the general momentum trend
878            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
879            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
880         END DO
881      END DO
882
883      ! ----------------
884      !  interior value   (2=<jk=<jpkm1)
885      ! ----------------
886#if defined key_z_first
887      DO jj = 2, jpjm1 
888         DO ji = 2, jpim1
889            DO jk = 2, jpkm1
890#else
891      DO jk = 2, jpkm1
892         DO jj = 2, jpjm1 
893            DO ji = fs_2, fs_jpim1   ! vector opt.
894#endif
895               ! hydrostatic pressure gradient along s-surfaces
896               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
897                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
898                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) / e1u(ji,jj)
899               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
900                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
901                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj)
902               ! add to the general momentum trend
903               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
904               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
905            END DO
906         END DO
907      END DO
908      !
909      IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) )   &
910         CALL ctl_stop('dyn:hpg_djc: failed to release workspace arrays')
911      !
912   END SUBROUTINE hpg_djc
913
914
915   SUBROUTINE hpg_rot( kt )
916      !!---------------------------------------------------------------------
917      !!                  ***  ROUTINE hpg_rot  ***
918      !!
919      !! ** Method  :   rotated axes scheme (Thiem and Berntsen 2005)
920      !!
921      !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005.
922      !!----------------------------------------------------------------------
923      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
924      USE oce     , ONLY:   zhpi    => ta       , zhpj    => sa       ! (ta,sa) used as 3D workspace
925      USE wrk_nemo, ONLY:   zdistr  => wrk_2d_1 , zsina   => wrk_2d_2 , zcosa  => wrk_2d_3
926      USE wrk_nemo, ONLY:   zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2
927      USE wrk_nemo, ONLY:   zhpitra => wrk_3d_3 , zhpine  => wrk_3d_4
928      USE wrk_nemo, ONLY:   zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6
929      USE wrk_nemo, ONLY:   zhpjtra => wrk_3d_7 , zhpjne  => wrk_3d_8
930      !! DCSE_NEMO: need additional directives for renamed module variables
931!FTRANS zhpi :I :I :z
932!FTRANS zhpj :I :I :z
933!FTRANS zhpiorg :I :I :z
934!FTRANS zhpirot :I :I :z
935!FTRANS zhpitra :I :I :z
936!FTRANS zhpine :I :I :z
937!FTRANS zhpjorg :I :I :z
938!FTRANS zhpjrot :I :I :z
939!FTRANS zhpjtra :I :I :z
940!FTRANS zhpjne :I :I :z
941      !!
942      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
943      !!
944      INTEGER  ::   ji, jj, jk          ! dummy loop indices
945      REAL(wp) ::   zforg, zcoef0, zuap, zmskd1, zmskd1m   ! temporary scalar
946      REAL(wp) ::   zfrot        , zvap, zmskd2, zmskd2m   !    "         "
947      !!----------------------------------------------------------------------
948
949      IF( wrk_in_use(2, 1,2,3)             .OR.   &
950          wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN
951         CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable')   ;   RETURN
952      ENDIF
953
954      IF( kt == nit000 ) THEN
955         IF(lwp) WRITE(numout,*)
956         IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend'
957         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, rotated axes scheme used'
958      ENDIF
959
960      ! -------------------------------
961      !  Local constant initialization
962      ! -------------------------------
963      zcoef0 = - grav * 0.5_wp
964      zforg  = 0.95_wp
965      zfrot  = 1._wp - zforg
966
967      ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance)
968      zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) )
969
970      ! sinus and cosinus of diagonal angle at F-point
971      zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) )
972      zcosa(:,:) = COS( zsina(:,:) )
973      zsina(:,:) = SIN( zsina(:,:) )
974
975      ! ---------------
976      !  Surface value
977      ! ---------------
978      ! compute and add to the general trend the pressure gradients along the axes
979      DO jj = 2, jpjm1
980         DO ji = fs_2, fs_jpim1   ! vector opt.
981            ! hydrostatic pressure gradient along s-surfaces
982            zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,1) * rhd(ji+1,jj,1)   &
983               &                                      - fse3t(ji  ,jj,1) * rhd(ji  ,jj,1)  )
984            zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,1) * rhd(ji,jj+1,1)   &
985               &                                      - fse3t(ji,jj  ,1) * rhd(ji,jj  ,1)  )
986            ! s-coordinate pressure gradient correction
987            zuap = -zcoef0 * ( rhd   (ji+1,jj  ,1) + rhd   (ji,jj,1) )   &
988               &           * ( fsdept(ji+1,jj  ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
989            zvap = -zcoef0 * ( rhd   (ji  ,jj+1,1) + rhd   (ji,jj,1) )   &
990               &           * ( fsdept(ji  ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
991            ! add to the general momentum trend
992            ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap )
993            va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap )
994         END DO
995      END DO
996
997      ! compute the pressure gradients in the diagonal directions
998      DO jj = 1, jpjm1
999         DO ji = 1, fs_jpim1   ! vector opt.
1000#if defined key_z_first
1001            zmskd1 = tmask_1(ji+1,jj+1) * tmask_1(ji  ,jj)      ! mask in the 1st diagnonal
1002            zmskd2 = tmask_1(ji  ,jj+1) * tmask_1(ji+1,jj)      ! mask in the 2nd diagnonal
1003#else
1004            zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji  ,jj,1)      ! mask in the 1st diagnonal
1005            zmskd2 = tmask(ji  ,jj+1,1) * tmask(ji+1,jj,1)      ! mask in the 2nd diagnonal
1006#endif
1007            ! hydrostatic pressure gradient along s-surfaces
1008            zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * (  fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1)   &
1009               &                                         - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  )
1010            zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * (  fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   &
1011               &                                         - fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  )
1012            ! s-coordinate pressure gradient correction
1013            zuap = -zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,1) + rhd   (ji  ,jj,1) )   &
1014               &                           * ( fsdept(ji+1,jj+1,1) - fsdept(ji  ,jj,1) )
1015            zvap = -zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,1) + rhd   (ji+1,jj,1) )   &
1016               &                           * ( fsdept(ji  ,jj+1,1) - fsdept(ji+1,jj,1) )
1017            ! back rotation
1018            zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   &
1019               &            - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
1020            zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   &
1021               &            + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
1022         END DO
1023      END DO
1024
1025      ! interpolate and add to the general trend the diagonal gradient
1026      DO jj = 2, jpjm1
1027         DO ji = fs_2, fs_jpim1   ! vector opt.
1028            ! averaging
1029            zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji  ,jj-1,1) )
1030            zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj  ,1) )
1031            ! add to the general momentum trend
1032            ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1) 
1033            va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1) 
1034         END DO
1035      END DO
1036
1037      ! -----------------
1038      ! 2. interior value (2=<jk=<jpkm1)
1039      ! -----------------
1040      ! compute and add to the general trend the pressure gradients along the axes
1041#if defined key_z_first
1042      DO jj = 2, jpjm1
1043         DO ji = 2, jpim1
1044            DO jk = 2, jpkm1
1045#else
1046      DO jk = 2, jpkm1
1047         DO jj = 2, jpjm1
1048            DO ji = fs_2, fs_jpim1   ! vector opt.
1049#endif
1050               ! hydrostatic pressure gradient along s-surfaces
1051               zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1)                                                 &
1052                  &              +  zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk  )   &
1053                  &                                        - fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk  )   &
1054                  &                                        + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   &
1055                  &                                        - fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1)  )
1056               zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1)                                                 &
1057                  &              +  zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk  )   &
1058                  &                                        - fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk  )   &
1059                  &                                        + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   &
1060                  &                                        - fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1)  )
1061               ! s-coordinate pressure gradient correction
1062               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   &
1063                  &            * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
1064               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   &
1065                  &            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
1066               ! add to the general momentum trend
1067               ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap )
1068               va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap )
1069            END DO
1070         END DO
1071      END DO
1072
1073      ! compute the pressure gradients in the diagonal directions
1074#if defined key_z_first
1075      DO jj = 1, jpjm1
1076         DO ji = 1, jpim1
1077            DO jk = 2, jpkm1
1078#else
1079      DO jk = 2, jpkm1
1080         DO jj = 1, jpjm1
1081            DO ji = 1, fs_jpim1   ! vector opt.
1082#endif
1083               zmskd1  = tmask(ji+1,jj+1,jk  ) * tmask(ji  ,jj,jk  )      ! level jk   mask in the 1st diagnonal
1084               zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji  ,jj,jk-1)      ! level jk-1    "               "     
1085               zmskd2  = tmask(ji  ,jj+1,jk  ) * tmask(ji+1,jj,jk  )      ! level jk   mask in the 2nd diagnonal
1086               zmskd2m = tmask(ji  ,jj+1,jk-1) * tmask(ji+1,jj,jk-1)      ! level jk-1    "               "     
1087               ! hydrostatic pressure gradient along s-surfaces
1088               zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1)                                                       &
1089                  &              + zdistr(ji,jj) * zmskd1  * ( fse3t(ji+1,jj+1,jk  ) * rhd(ji+1,jj+1,jk)     &
1090                  &                                           -fse3t(ji  ,jj  ,jk  ) * rhd(ji  ,jj  ,jk) )   &
1091                  &              + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1)   &
1092                  &                                           -fse3t(ji  ,jj  ,jk-1) * rhd(ji  ,jj  ,jk-1) )
1093               zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1)                                                       &
1094                  &              + zdistr(ji,jj) * zmskd2  * ( fse3t(ji  ,jj+1,jk  ) * rhd(ji  ,jj+1,jk)     &
1095                  &                                           -fse3t(ji+1,jj  ,jk  ) * rhd(ji+1,jj,  jk) )   &
1096                  &              + zdistr(ji,jj) * zmskd2m * ( fse3t(ji  ,jj+1,jk-1) * rhd(ji  ,jj+1,jk-1)   &
1097                  &                                           -fse3t(ji+1,jj  ,jk-1) * rhd(ji+1,jj,  jk-1) )
1098               ! s-coordinate pressure gradient correction
1099               zuap = - zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,jk) + rhd   (ji  ,jj,jk) )   &
1100                  &                            * ( fsdept(ji+1,jj+1,jk) - fsdept(ji  ,jj,jk) )
1101               zvap = - zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji+1,jj,jk) )   &
1102                  &                            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji+1,jj,jk) )
1103               ! back rotation
1104               zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   &
1105                  &             - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
1106               zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   &
1107                  &             + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
1108            END DO
1109         END DO
1110      END DO
1111
1112      ! interpolate and add to the general trend
1113#if defined key_z_first
1114      DO jj = 2, jpjm1
1115         DO ji = 2, jpim1
1116            DO jk = 2, jpkm1
1117#else
1118      DO jk = 2, jpkm1
1119         DO jj = 2, jpjm1
1120            DO ji = fs_2, fs_jpim1   ! vector opt.
1121#endif
1122               ! averaging
1123               zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji  ,jj-1,jk) )
1124               zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj  ,jk) )
1125               ! add to the general momentum trend
1126               ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk) 
1127               va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk) 
1128            END DO
1129         END DO
1130      END DO
1131      !
1132      IF( wrk_not_released(2, 1,2,3)           .OR.   &
1133          wrk_not_released(3, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays')
1134      !
1135   END SUBROUTINE hpg_rot
1136
1137   !!======================================================================
1138END MODULE dynhpg
Note: See TracBrowser for help on using the repository browser.