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 trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 503

Last change on this file since 503 was 503, checked in by opalod, 18 years ago

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