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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • Property svn:keywords set to Id
File size: 52.2 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   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
17   !!                 !           (A. Coward) suppression of hel, wdj and rot options
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dyn_hpg      : update the momentum trend with the now horizontal
22   !!                  gradient of the hydrostatic pressure
23   !!   dyn_hpg_init : initialisation and control of options
24   !!       hpg_zco  : z-coordinate scheme
25   !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
26   !!       hpg_sco  : s-coordinate (standard jacobian formulation)
27   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
28   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial)
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   USE wrk_nemo        ! Memory Allocation
40   USE timing          ! Timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   dyn_hpg        ! routine called by step module
46   PUBLIC   dyn_hpg_init   ! routine called by opa module
47
48   !                                    !!* Namelist namdyn_hpg : hydrostatic pressure gradient
49   LOGICAL , PUBLIC ::   ln_hpg_zco      !: z-coordinate - full steps
50   LOGICAL , PUBLIC ::   ln_hpg_zps      !: z-coordinate - partial steps (interpolation)
51   LOGICAL , PUBLIC ::   ln_hpg_sco      !: s-coordinate (standard jacobian formulation)
52   LOGICAL , PUBLIC ::   ln_hpg_djc      !: s-coordinate (Density Jacobian with Cubic polynomial)
53   LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme)
54   LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag
55
56   INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE dyn_hpg( kt )
69      !!---------------------------------------------------------------------
70      !!                  ***  ROUTINE dyn_hpg  ***
71      !!
72      !! ** Method  :   Call the hydrostatic pressure gradient routine
73      !!              using the scheme defined in the namelist
74      !!
75      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
76      !!             - Save the trend (l_trddyn=T)
77      !!----------------------------------------------------------------------
78      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
79      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
80      !!----------------------------------------------------------------------
81      !
82      IF( nn_timing == 1 )  CALL timing_start('dyn_hpg')
83      !
84      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
85         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
86         ztrdu(:,:,:) = ua(:,:,:)
87         ztrdv(:,:,:) = va(:,:,:)
88      ENDIF
89      !
90      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
91      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
92      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
93      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
94      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
95      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
96      END SELECT
97      !
98      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
99         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
100         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
101         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
102         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
103      ENDIF
104      !
105      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
106         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
107      !
108      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg')
109      !
110   END SUBROUTINE dyn_hpg
111
112
113   SUBROUTINE dyn_hpg_init
114      !!----------------------------------------------------------------------
115      !!                 ***  ROUTINE dyn_hpg_init  ***
116      !!
117      !! ** Purpose :   initializations for the hydrostatic pressure gradient
118      !!              computation and consistency control
119      !!
120      !! ** Action  :   Read the namelist namdyn_hpg and check the consistency
121      !!      with the type of vertical coordinate used (zco, zps, sco)
122      !!----------------------------------------------------------------------
123      INTEGER ::   ioptio = 0      ! temporary integer
124      INTEGER ::   ios             ! Local integer output status for namelist read
125      !!
126      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
127         &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp
128      !!----------------------------------------------------------------------
129      !
130      REWIND( numnam_ref )              ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient
131      READ  ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901)
132901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )
133
134      REWIND( numnam_cfg )              ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient
135      READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
136902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )
137      IF(lwm) WRITE ( numond, namdyn_hpg )
138      !
139      IF(lwp) THEN                   ! Control print
140         WRITE(numout,*)
141         WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
142         WRITE(numout,*) '~~~~~~~~~~~~'
143         WRITE(numout,*) '   Namelist namdyn_hpg : choice of hpg scheme'
144         WRITE(numout,*) '      z-coord. - full steps                             ln_hpg_zco    = ', ln_hpg_zco
145         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps
146         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco
147         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
148         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj
149         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp
150      ENDIF
151      !
152      IF( ln_hpg_djc )   &
153         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method &
154                           & currently disabled (bugs under investigation). Please select &
155                           & either  ln_hpg_sco or  ln_hpg_prj instead')
156      !
157      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   &
158         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&
159                           & the standard jacobian formulation hpg_sco or &
160                           & the pressure jacobian formulation hpg_prj')
161      !
162      !                               ! Set nhpg from ln_hpg_... flags
163      IF( ln_hpg_zco )   nhpg = 0
164      IF( ln_hpg_zps )   nhpg = 1
165      IF( ln_hpg_sco )   nhpg = 2
166      IF( ln_hpg_djc )   nhpg = 3
167      IF( ln_hpg_prj )   nhpg = 4
168      !
169      !                               ! Consistency check
170      ioptio = 0
171      IF( ln_hpg_zco )   ioptio = ioptio + 1
172      IF( ln_hpg_zps )   ioptio = ioptio + 1
173      IF( ln_hpg_sco )   ioptio = ioptio + 1
174      IF( ln_hpg_djc )   ioptio = ioptio + 1
175      IF( ln_hpg_prj )   ioptio = ioptio + 1
176      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
177      !
178   END SUBROUTINE dyn_hpg_init
179
180
181   SUBROUTINE hpg_zco( kt )
182      !!---------------------------------------------------------------------
183      !!                  ***  ROUTINE hpg_zco  ***
184      !!
185      !! ** Method  :   z-coordinate case, levels are horizontal surfaces.
186      !!      The now hydrostatic pressure gradient at a given level, jk,
187      !!      is computed by taking the vertical integral of the in-situ
188      !!      density gradient along the model level from the suface to that
189      !!      level:    zhpi = grav .....
190      !!                zhpj = grav .....
191      !!      add it to the general momentum trend (ua,va).
192      !!            ua = ua - 1/e1u * zhpi
193      !!            va = va - 1/e2v * zhpj
194      !!
195      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
196      !!----------------------------------------------------------------------
197      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
198      !!
199      INTEGER  ::   ji, jj, jk       ! dummy loop indices
200      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
201      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
202      !!----------------------------------------------------------------------
203      !
204      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
205      !
206      IF( kt == nit000 ) THEN
207         IF(lwp) WRITE(numout,*)
208         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
209         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
210      ENDIF
211
212      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
213
214      ! Surface value
215      DO jj = 2, jpjm1
216         DO ji = fs_2, fs_jpim1   ! vector opt.
217            zcoef1 = zcoef0 * fse3w(ji,jj,1)
218            ! hydrostatic pressure gradient
219            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
220            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
221            ! add to the general momentum trend
222            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
223            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
224         END DO
225      END DO
226
227      !
228      ! interior value (2=<jk=<jpkm1)
229      DO jk = 2, jpkm1
230         DO jj = 2, jpjm1
231            DO ji = fs_2, fs_jpim1   ! vector opt.
232               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
233               ! hydrostatic pressure gradient
234               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
235                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
236                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
237
238               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
239                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
240                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
241               ! add to the general momentum trend
242               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
243               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
244            END DO
245         END DO
246      END DO
247      !
248      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
249      !
250   END SUBROUTINE hpg_zco
251
252
253   SUBROUTINE hpg_zps( kt )
254      !!---------------------------------------------------------------------
255      !!                 ***  ROUTINE hpg_zps  ***
256      !!
257      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
258      !!
259      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
260      !!----------------------------------------------------------------------
261      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
262      !!
263      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
264      INTEGER  ::   iku, ikv                         ! temporary integers
265      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
266      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
267      !!----------------------------------------------------------------------
268      !
269      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
270      !
271      IF( kt == nit000 ) THEN
272         IF(lwp) WRITE(numout,*)
273         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
274         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
275      ENDIF
276
277
278      ! Local constant initialization
279      zcoef0 = - grav * 0.5_wp
280
281      !  Surface value (also valid in partial step case)
282      DO jj = 2, jpjm1
283         DO ji = fs_2, fs_jpim1   ! vector opt.
284            zcoef1 = zcoef0 * fse3w(ji,jj,1)
285            ! hydrostatic pressure gradient
286            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
287            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
288            ! add to the general momentum trend
289            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
290            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
291         END DO
292      END DO
293
294
295      ! interior value (2=<jk=<jpkm1)
296      DO jk = 2, jpkm1
297         DO jj = 2, jpjm1
298            DO ji = fs_2, fs_jpim1   ! vector opt.
299               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
300               ! hydrostatic pressure gradient
301               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
302                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
303                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
304
305               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
306                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
307                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
308               ! add to the general momentum trend
309               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
310               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
311            END DO
312         END DO
313      END DO
314
315
316      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
317# if defined key_vectopt_loop
318         jj = 1
319         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
320# else
321      DO jj = 2, jpjm1
322         DO ji = 2, jpim1
323# endif
324            iku = mbku(ji,jj)
325            ikv = mbkv(ji,jj)
326            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) )
327            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) )
328            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
329               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
330               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
331                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
332               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
333            ENDIF
334            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
335               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
336               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
337                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
338               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
339            ENDIF
340# if ! defined key_vectopt_loop
341         END DO
342# endif
343      END DO
344      !
345      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
346      !
347   END SUBROUTINE hpg_zps
348
349
350   SUBROUTINE hpg_sco( kt )
351      !!---------------------------------------------------------------------
352      !!                  ***  ROUTINE hpg_sco  ***
353      !!
354      !! ** Method  :   s-coordinate case. Jacobian scheme.
355      !!      The now hydrostatic pressure gradient at a given level, jk,
356      !!      is computed by taking the vertical integral of the in-situ
357      !!      density gradient along the model level from the suface to that
358      !!      level. s-coordinates (ln_sco): a corrective term is added
359      !!      to the horizontal pressure gradient :
360      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
361      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
362      !!      add it to the general momentum trend (ua,va).
363      !!         ua = ua - 1/e1u * zhpi
364      !!         va = va - 1/e2v * zhpj
365      !!
366      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
367      !!----------------------------------------------------------------------
368      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
369      !!
370      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
371      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
372      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
373      !!----------------------------------------------------------------------
374      !
375      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
376      !
377      IF( kt == nit000 ) THEN
378         IF(lwp) WRITE(numout,*)
379         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
380         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
381      ENDIF
382
383      ! Local constant initialization
384      zcoef0 = - grav * 0.5_wp
385      ! To use density and not density anomaly
386      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
387      ELSE                 ;     znad = 0._wp         ! Fixed volume
388      ENDIF
389
390      ! Surface value
391      DO jj = 2, jpjm1
392         DO ji = fs_2, fs_jpim1   ! vector opt.
393            ! hydrostatic pressure gradient along s-surfaces
394            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
395               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
396            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
397               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
398            ! s-coordinate pressure gradient correction
399            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
400               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
401            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
402               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
403            ! add to the general momentum trend
404            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
405            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
406         END DO
407      END DO
408
409      ! interior value (2=<jk=<jpkm1)
410      DO jk = 2, jpkm1
411         DO jj = 2, jpjm1
412            DO ji = fs_2, fs_jpim1   ! vector opt.
413               ! hydrostatic pressure gradient along s-surfaces
414               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
415                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
416                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
417               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
418                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
419                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
420               ! s-coordinate pressure gradient correction
421               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
422                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
423               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
424                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
425               ! add to the general momentum trend
426               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
427               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
428            END DO
429         END DO
430      END DO
431      !
432      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
433      !
434   END SUBROUTINE hpg_sco
435
436   SUBROUTINE hpg_djc( kt )
437      !!---------------------------------------------------------------------
438      !!                  ***  ROUTINE hpg_djc  ***
439      !!
440      !! ** Method  :   Density Jacobian with Cubic polynomial scheme
441      !!
442      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
443      !!----------------------------------------------------------------------
444      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
445      !!
446      INTEGER  ::   ji, jj, jk          ! dummy loop indices
447      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
448      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
449      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
450      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
451      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw
452      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow
453      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k
454      !!----------------------------------------------------------------------
455      !
456      CALL wrk_alloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )
457      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
458      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
459      !
460
461      IF( kt == nit000 ) THEN
462         IF(lwp) WRITE(numout,*)
463         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
464         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
465      ENDIF
466
467      ! Local constant initialization
468      zcoef0 = - grav * 0.5_wp
469      z1_10  = 1._wp / 10._wp
470      z1_12  = 1._wp / 12._wp
471
472      !----------------------------------------------------------------------------------------
473      !  compute and store in provisional arrays elementary vertical and horizontal differences
474      !----------------------------------------------------------------------------------------
475
476!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
477
478      DO jk = 2, jpkm1
479         DO jj = 2, jpjm1
480            DO ji = fs_2, fs_jpim1   ! vector opt.
481               drhoz(ji,jj,jk) = rhd   (ji  ,jj  ,jk) - rhd   (ji,jj,jk-1)
482               dzz  (ji,jj,jk) = fsde3w(ji  ,jj  ,jk) - fsde3w(ji,jj,jk-1)
483               drhox(ji,jj,jk) = rhd   (ji+1,jj  ,jk) - rhd   (ji,jj,jk  )
484               dzx  (ji,jj,jk) = fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk  )
485               drhoy(ji,jj,jk) = rhd   (ji  ,jj+1,jk) - rhd   (ji,jj,jk  )
486               dzy  (ji,jj,jk) = fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk  )
487            END DO
488         END DO
489      END DO
490
491      !-------------------------------------------------------------------------
492      ! compute harmonic averages using eq. 5.18
493      !-------------------------------------------------------------------------
494      zep = 1.e-15
495
496!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
497!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
498
499      DO jk = 2, jpkm1
500         DO jj = 2, jpjm1
501            DO ji = fs_2, fs_jpim1   ! vector opt.
502               cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
503
504               cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  )
505               cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  )
506
507               cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  )
508               cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  )
509
510               IF( cffw > zep) THEN
511                  drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   &
512                     &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
513               ELSE
514                  drhow(ji,jj,jk) = 0._wp
515               ENDIF
516
517               dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   &
518                  &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
519
520               IF( cffu > zep ) THEN
521                  drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   &
522                     &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
523               ELSE
524                  drhou(ji,jj,jk ) = 0._wp
525               ENDIF
526
527               IF( cffx > zep ) THEN
528                  dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   &
529                     &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
530               ELSE
531                  dzu(ji,jj,jk) = 0._wp
532               ENDIF
533
534               IF( cffv > zep ) THEN
535                  drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   &
536                     &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
537               ELSE
538                  drhov(ji,jj,jk) = 0._wp
539               ENDIF
540
541               IF( cffy > zep ) THEN
542                  dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   &
543                     &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
544               ELSE
545                  dzv(ji,jj,jk) = 0._wp
546               ENDIF
547
548            END DO
549         END DO
550      END DO
551
552      !----------------------------------------------------------------------------------
553      ! apply boundary conditions at top and bottom using 5.36-5.37
554      !----------------------------------------------------------------------------------
555      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  )
556      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  )
557      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  )
558
559      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
560      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
561      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
562
563
564      !--------------------------------------------------------------
565      ! Upper half of top-most grid box, compute and store
566      !-------------------------------------------------------------
567
568!!bug gm   :  e3w-de3w = 0.5*e3w  ....  and de3w(2)-de3w(1)=e3w(2) ....   to be verified
569!          true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
570
571      DO jj = 2, jpjm1
572         DO ji = fs_2, fs_jpim1   ! vector opt.
573            rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )               &
574               &                   * (  rhd(ji,jj,1)                                    &
575               &                     + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         &
576               &                              * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   &
577               &                              / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  )
578         END DO
579      END DO
580
581!!bug gm    : here also, simplification is possible
582!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
583
584      DO jk = 2, jpkm1
585         DO jj = 2, jpjm1
586            DO ji = fs_2, fs_jpim1   ! vector opt.
587
588               rho_k(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj,jk) + rhd   (ji,jj,jk-1) )                                   &
589                  &                     * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) )                                   &
590                  &            - grav * z1_10 * (                                                                     &
591                  &     ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) )                                                     &
592                  &   * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
593                  &   - ( dzw   (ji,jj,jk) - dzw   (ji,jj,jk-1) )                                                     &
594                  &   * ( rhd   (ji,jj,jk) - rhd   (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   &
595                  &                             )
596
597               rho_i(ji,jj,jk) = zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )                                   &
598                  &                     * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) )                                   &
599                  &            - grav* z1_10 * (                                                                      &
600                  &     ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) )                                                     &
601                  &   * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
602                  &   - ( dzu   (ji+1,jj,jk) - dzu   (ji,jj,jk) )                                                     &
603                  &   * ( rhd   (ji+1,jj,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   &
604                  &                            )
605
606               rho_j(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )                                   &
607                  &                     * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) )                                   &
608                  &            - grav* z1_10 * (                                                                      &
609                  &     ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) )                                                     &
610                  &   * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
611                  &   - ( dzv   (ji,jj+1,jk) - dzv   (ji,jj,jk) )                                                     &
612                  &   * ( rhd   (ji,jj+1,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   &
613                  &                            )
614
615            END DO
616         END DO
617      END DO
618      CALL lbc_lnk(rho_k,'W',1.)
619      CALL lbc_lnk(rho_i,'U',1.)
620      CALL lbc_lnk(rho_j,'V',1.)
621
622
623      ! ---------------
624      !  Surface value
625      ! ---------------
626      DO jj = 2, jpjm1
627         DO ji = fs_2, fs_jpim1   ! vector opt.
628            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
629            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
630            ! add to the general momentum trend
631            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
632            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
633         END DO
634      END DO
635
636      ! ----------------
637      !  interior value   (2=<jk=<jpkm1)
638      ! ----------------
639      DO jk = 2, jpkm1
640         DO jj = 2, jpjm1
641            DO ji = fs_2, fs_jpim1   ! vector opt.
642               ! hydrostatic pressure gradient along s-surfaces
643               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
644                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
645                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) / e1u(ji,jj)
646               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
647                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
648                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj)
649               ! add to the general momentum trend
650               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
651               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
652            END DO
653         END DO
654      END DO
655      !
656      CALL wrk_dealloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )
657      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
658      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
659      !
660   END SUBROUTINE hpg_djc
661
662
663   SUBROUTINE hpg_prj( kt )
664      !!---------------------------------------------------------------------
665      !!                  ***  ROUTINE hpg_prj  ***
666      !!
667      !! ** Method  :   s-coordinate case.
668      !!      A Pressure-Jacobian horizontal pressure gradient method
669      !!      based on the constrained cubic-spline interpolation for
670      !!      all vertical coordinate systems
671      !!
672      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
673      !!             - Save the trend (l_trddyn=T)
674      !!
675      !!----------------------------------------------------------------------
676      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
677      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
678      !!
679      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
680      REAL(wp) ::   zcoef0, znad                    ! temporary scalars
681      !!
682      !! The local variables for the correction term
683      INTEGER  :: jk1, jis, jid, jjs, jjd
684      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
685      REAL(wp) :: zrhdt1
686      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
687      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
688      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
689      !!----------------------------------------------------------------------
690      !
691      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
692      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
693      !
694      IF( kt == nit000 ) THEN
695         IF(lwp) WRITE(numout,*)
696         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
697         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
698      ENDIF
699
700      !!----------------------------------------------------------------------
701      ! Local constant initialization
702      zcoef0 = - grav
703      znad = 0.0_wp
704      IF( lk_vvl ) znad = 1._wp
705
706      ! Clean 3-D work arrays
707      zhpi(:,:,:) = 0._wp
708      zrhh(:,:,:) = rhd(:,:,:)
709
710      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
711      DO jj = 1, jpj
712        DO ji = 1, jpi
713          jk = mbathy(ji,jj)
714          IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp
715          ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
716          ELSE IF(jk < jpkm1) THEN
717             DO jkk = jk+1, jpk
718                zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), &
719                                         fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
720             END DO
721          ENDIF
722        END DO
723      END DO
724
725      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
726      DO jj = 1, jpj;   DO ji = 1, jpi
727          zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad
728      END DO        ;   END DO
729
730      DO jk = 2, jpk;   DO jj = 1, jpj;   DO ji = 1, jpi
731          zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)
732      END DO        ;   END DO        ;   END DO
733
734      fsp(:,:,:) = zrhh(:,:,:)
735      xsp(:,:,:) = zdept(:,:,:)
736
737      ! Construct the vertical density profile with the
738      ! constrained cubic spline interpolation
739      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3
740      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)
741
742      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
743      DO jj = 2, jpj
744        DO ji = 2, jpi
745          zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), &
746                                         bsp(ji,jj,1),   csp(ji,jj,1), &
747                                         dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1)
748
749          ! assuming linear profile across the top half surface layer
750          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * zrhdt1
751        END DO
752      END DO
753
754      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
755      DO jk = 2, jpkm1
756        DO jj = 2, jpj
757          DO ji = 2, jpi
758            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          &
759                             integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&
760                                    asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), &
761                                    csp(ji,jj,jk-1),    dsp(ji,jj,jk-1))
762          END DO
763        END DO
764      END DO
765
766      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
767      DO jj = 2, jpjm1
768        DO ji = 2, jpim1
769          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation
770          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation
771        END DO
772      END DO
773
774      DO jk = 2, jpkm1
775        DO jj = 2, jpjm1
776          DO ji = 2, jpim1
777            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk)
778            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk)
779          END DO
780        END DO
781      END DO
782
783      DO jk = 1, jpkm1
784        DO jj = 2, jpjm1
785          DO ji = 2, jpim1
786            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)
787            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)
788          END DO
789        END DO
790      END DO
791
792      DO jk = 1, jpkm1
793        DO jj = 2, jpjm1
794          DO ji = 2, jpim1
795            zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
796            zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
797            zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
798            zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
799          END DO
800        END DO
801      END DO
802
803
804      DO jk = 1, jpkm1
805        DO jj = 2, jpjm1
806          DO ji = 2, jpim1
807            zpwes = 0._wp; zpwed = 0._wp
808            zpnss = 0._wp; zpnsd = 0._wp
809            zuijk = zu(ji,jj,jk)
810            zvijk = zv(ji,jj,jk)
811
812            !!!!!     for u equation
813            IF( jk <= mbku(ji,jj) ) THEN
814               IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
815                 jis = ji + 1; jid = ji
816               ELSE
817                 jis = ji;     jid = ji +1
818               ENDIF
819
820               ! integrate the pressure on the shallow side
821               jk1 = jk
822               DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
823                 IF( jk1 == mbku(ji,jj) ) THEN
824                   zuijk = -zdept(jis,jj,jk1)
825                   EXIT
826                 ENDIF
827                 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
828                 zpwes = zpwes +                                    &
829                      integ_spline(zdept(jis,jj,jk1), zdeps,            &
830                             asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
831                             csp(jis,jj,jk1),    dsp(jis,jj,jk1))
832                 jk1 = jk1 + 1
833               END DO
834
835               ! integrate the pressure on the deep side
836               jk1 = jk
837               DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
838                 IF( jk1 == 1 ) THEN
839                   zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
840                   zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
841                                                     bsp(jid,jj,1),   csp(jid,jj,1), &
842                                                     dsp(jid,jj,1)) * zdeps
843                   zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
844                   EXIT
845                 ENDIF
846                 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
847                 zpwed = zpwed +                                        &
848                        integ_spline(zdeps,              zdept(jid,jj,jk1), &
849                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
850                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
851                 jk1 = jk1 - 1
852               END DO
853
854               ! update the momentum trends in u direction
855
856               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
857               IF( lk_vvl ) THEN
858                 zdpdx2 = zcoef0 / e1u(ji,jj) * &
859                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
860                ELSE
861                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
862               ENDIF
863
864               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * &
865               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
866            ENDIF
867
868            !!!!!     for v equation
869            IF( jk <= mbkv(ji,jj) ) THEN
870               IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
871                 jjs = jj + 1; jjd = jj
872               ELSE
873                 jjs = jj    ; jjd = jj + 1
874               ENDIF
875
876               ! integrate the pressure on the shallow side
877               jk1 = jk
878               DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
879                 IF( jk1 == mbkv(ji,jj) ) THEN
880                   zvijk = -zdept(ji,jjs,jk1)
881                   EXIT
882                 ENDIF
883                 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
884                 zpnss = zpnss +                                      &
885                        integ_spline(zdept(ji,jjs,jk1), zdeps,            &
886                               asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
887                               csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
888                 jk1 = jk1 + 1
889               END DO
890
891               ! integrate the pressure on the deep side
892               jk1 = jk
893               DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
894                 IF( jk1 == 1 ) THEN
895                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
896                   zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
897                                                     bsp(ji,jjd,1),   csp(ji,jjd,1), &
898                                                     dsp(ji,jjd,1) ) * zdeps
899                   zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
900                   EXIT
901                 ENDIF
902                 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
903                 zpnsd = zpnsd +                                        &
904                        integ_spline(zdeps,              zdept(ji,jjd,jk1), &
905                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
906                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
907                 jk1 = jk1 - 1
908               END DO
909
910
911               ! update the momentum trends in v direction
912
913               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
914               IF( lk_vvl ) THEN
915                   zdpdy2 = zcoef0 / e2v(ji,jj) * &
916                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
917               ELSE
918                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
919               ENDIF
920
921               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*&
922               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
923            ENDIF
924
925
926           END DO
927        END DO
928      END DO
929      !
930      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
931      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
932      !
933   END SUBROUTINE hpg_prj
934
935   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
936      !!----------------------------------------------------------------------
937      !!                 ***  ROUTINE cspline  ***
938      !!
939      !! ** Purpose :   constrained cubic spline interpolation
940      !!
941      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3
942      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
943      !!
944      !!----------------------------------------------------------------------
945      IMPLICIT NONE
946      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate
947      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of
948                                                                    ! the interpoated function
949      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline
950                                                                    ! 2: Linear
951
952      ! Local Variables
953      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
954      INTEGER  ::   jpi, jpj, jpkm1
955      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
956      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha
957      REAL(wp) ::   zdf(size(fsp,3))
958      !!----------------------------------------------------------------------
959
960      jpi   = size(fsp,1)
961      jpj   = size(fsp,2)
962      jpkm1 = size(fsp,3) - 1
963
964
965      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
966         DO ji = 1, jpi
967            DO jj = 1, jpj
968           !!Fritsch&Butland's method, 1984 (preferred, but more computation)
969           !    DO jk = 2, jpkm1-1
970           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)
971           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
972           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
973           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
974           !
975           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
976           !
977           !       IF(zdf1 * zdf2 <= 0._wp) THEN
978           !           zdf(jk) = 0._wp
979           !       ELSE
980           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
981           !       ENDIF
982           !    END DO
983
984           !!Simply geometric average
985               DO jk = 2, jpkm1-1
986                  zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))
987                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))
988
989                  IF(zdf1 * zdf2 <= 0._wp) THEN
990                     zdf(jk) = 0._wp
991                  ELSE
992                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
993                  ENDIF
994               END DO
995
996               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
997                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2)
998               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
999                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
1000                          & 0.5_wp * zdf(jpkm1 - 1)
1001
1002               DO jk = 1, jpkm1 - 1
1003                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1004                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1005                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1006                 zddf1  = -2._wp * ztmp1 + ztmp2
1007                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1008                 zddf2  =  2._wp * ztmp1 - ztmp2
1009
1010                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1011                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1012                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1013                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1014                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
1015                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
1016                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
1017                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
1018                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
1019               END DO
1020            END DO
1021         END DO
1022
1023      ELSE IF (polynomial_type == 2) THEN     ! Linear
1024         DO ji = 1, jpi
1025            DO jj = 1, jpj
1026               DO jk = 1, jpkm1-1
1027                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1028                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1029
1030                  dsp(ji,jj,jk) = 0._wp
1031                  csp(ji,jj,jk) = 0._wp
1032                  bsp(ji,jj,jk) = ztmp1 / zdxtmp
1033                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1034               END DO
1035            END DO
1036         END DO
1037
1038      ELSE
1039           CALL ctl_stop( 'invalid polynomial type in cspline' )
1040      ENDIF
1041
1042
1043   END SUBROUTINE cspline
1044
1045
1046   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)
1047      !!----------------------------------------------------------------------
1048      !!                 ***  ROUTINE interp1  ***
1049      !!
1050      !! ** Purpose :   1-d linear interpolation
1051      !!
1052      !! ** Method  :
1053      !!                interpolation is straight forward
1054      !!                extrapolation is also permitted (no value limit)
1055      !!
1056      !!----------------------------------------------------------------------
1057      IMPLICIT NONE
1058      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr
1059      REAL(wp)             ::  f ! result of the interpolation (extrapolation)
1060      REAL(wp)             ::  zdeltx
1061      !!----------------------------------------------------------------------
1062
1063      zdeltx = xr - xl
1064      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
1065        f = 0.5_wp * (fl + fr)
1066      ELSE
1067        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1068      ENDIF
1069
1070   END FUNCTION interp1
1071
1072   FUNCTION interp2(x, a, b, c, d)  RESULT(f)
1073      !!----------------------------------------------------------------------
1074      !!                 ***  ROUTINE interp1  ***
1075      !!
1076      !! ** Purpose :   1-d constrained cubic spline interpolation
1077      !!
1078      !! ** Method  :  cubic spline interpolation
1079      !!
1080      !!----------------------------------------------------------------------
1081      IMPLICIT NONE
1082      REAL(wp), INTENT(in) ::  x, a, b, c, d
1083      REAL(wp)             ::  f ! value from the interpolation
1084      !!----------------------------------------------------------------------
1085
1086      f = a + x* ( b + x * ( c + d * x ) )
1087
1088   END FUNCTION interp2
1089
1090
1091   FUNCTION interp3(x, a, b, c, d)  RESULT(f)
1092      !!----------------------------------------------------------------------
1093      !!                 ***  ROUTINE interp1  ***
1094      !!
1095      !! ** Purpose :   Calculate the first order of deriavtive of
1096      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3
1097      !!
1098      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2
1099      !!
1100      !!----------------------------------------------------------------------
1101      IMPLICIT NONE
1102      REAL(wp), INTENT(in) ::  x, a, b, c, d
1103      REAL(wp)             ::  f ! value from the interpolation
1104      !!----------------------------------------------------------------------
1105
1106      f = b + x * ( 2._wp * c + 3._wp * d * x)
1107
1108   END FUNCTION interp3
1109
1110
1111   FUNCTION integ_spline(xl, xr, a, b, c, d)  RESULT(f)
1112      !!----------------------------------------------------------------------
1113      !!                 ***  ROUTINE interp1  ***
1114      !!
1115      !! ** Purpose :   1-d constrained cubic spline integration
1116      !!
1117      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr
1118      !!
1119      !!----------------------------------------------------------------------
1120      IMPLICIT NONE
1121      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d
1122      REAL(wp)             ::  za1, za2, za3
1123      REAL(wp)             ::  f                   ! integration result
1124      !!----------------------------------------------------------------------
1125
1126      za1 = 0.5_wp * b
1127      za2 = c / 3.0_wp
1128      za3 = 0.25_wp * d
1129
1130      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1131         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1132
1133   END FUNCTION integ_spline
1134
1135
1136   !!======================================================================
1137END MODULE dynhpg
1138
Note: See TracBrowser for help on using the repository browser.