source: NEMO/trunk/src/OCE/DYN/dynhpg.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 8 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 69.0 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   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dyn_hpg      : update the momentum trend with the now horizontal
23   !!                  gradient of the hydrostatic pressure
24   !!   dyn_hpg_init : initialisation and control of options
25   !!       hpg_zco  : z-coordinate scheme
26   !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
27   !!       hpg_sco  : s-coordinate (standard jacobian formulation)
28   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf
29   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
30   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial)
31   !!----------------------------------------------------------------------
32   USE oce             ! ocean dynamics and tracers
33   USE isf_oce , ONLY : risfload  ! ice shelf  (risfload variable)
34   USE isfload , ONLY : isf_load  ! ice shelf  (isf_load routine )
35   USE sbc_oce         ! surface variable (only for the flag with ice shelf)
36   USE dom_oce         ! ocean space and time domain
37   USE wet_dry         ! wetting and drying
38   USE phycst          ! physical constants
39   USE trd_oce         ! trends: ocean variables
40   USE trddyn          ! trend manager: dynamics
41   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
42   !
43   USE in_out_manager  ! I/O manager
44   USE prtctl          ! Print control
45   USE lbclnk          ! lateral boundary condition
46   USE lib_mpp         ! MPP library
47   USE eosbn2          ! compute density
48   USE timing          ! Timing
49   USE iom
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   dyn_hpg        ! routine called by step module
55   PUBLIC   dyn_hpg_init   ! routine called by opa module
56
57   !                                !!* Namelist namdyn_hpg : hydrostatic pressure gradient
58   LOGICAL, PUBLIC ::   ln_hpg_zco   !: z-coordinate - full steps
59   LOGICAL, PUBLIC ::   ln_hpg_zps   !: z-coordinate - partial steps (interpolation)
60   LOGICAL, PUBLIC ::   ln_hpg_sco   !: s-coordinate (standard jacobian formulation)
61   LOGICAL, PUBLIC ::   ln_hpg_djc   !: s-coordinate (Density Jacobian with Cubic polynomial)
62   LOGICAL, PUBLIC ::   ln_hpg_prj   !: s-coordinate (Pressure Jacobian scheme)
63   LOGICAL, PUBLIC ::   ln_hpg_isf   !: s-coordinate similar to sco modify for isf
64
65   !                                !! Flag to control the type of hydrostatic pressure gradient
66   INTEGER, PARAMETER ::   np_ERROR  =-10   ! error in specification of lateral diffusion
67   INTEGER, PARAMETER ::   np_zco    =  0   ! z-coordinate - full steps
68   INTEGER, PARAMETER ::   np_zps    =  1   ! z-coordinate - partial steps (interpolation)
69   INTEGER, PARAMETER ::   np_sco    =  2   ! s-coordinate (standard jacobian formulation)
70   INTEGER, PARAMETER ::   np_djc    =  3   ! s-coordinate (Density Jacobian with Cubic polynomial)
71   INTEGER, PARAMETER ::   np_prj    =  4   ! s-coordinate (Pressure Jacobian scheme)
72   INTEGER, PARAMETER ::   np_isf    =  5   ! s-coordinate similar to sco modify for isf
73   !
74   INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
75
76   !! * Substitutions
77#  include "do_loop_substitute.h90"
78   !!----------------------------------------------------------------------
79   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
80   !! $Id$
81   !! Software governed by the CeCILL license (see ./LICENSE)
82   !!----------------------------------------------------------------------
83CONTAINS
84
85   SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs )
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 (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
93      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
94      !!----------------------------------------------------------------------
95      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
96      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
97      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
98      !
99      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
100      !!----------------------------------------------------------------------
101      !
102      IF( ln_timing )   CALL timing_start('dyn_hpg')
103      !
104      IF( l_trddyn ) THEN                    ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn)
105         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
106         ztrdu(:,:,:) = puu(:,:,:,Krhs)
107         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
108      ENDIF
109      !
110      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
111      CASE ( np_zco )   ;   CALL hpg_zco    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate
112      CASE ( np_zps )   ;   CALL hpg_zps    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate plus partial steps (interpolation)
113      CASE ( np_sco )   ;   CALL hpg_sco    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (standard jacobian formulation)
114      CASE ( np_djc )   ;   CALL hpg_djc    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Density Jacobian with Cubic polynomial)
115      CASE ( np_prj )   ;   CALL hpg_prj    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Pressure Jacobian scheme)
116      CASE ( np_isf )   ;   CALL hpg_isf    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate similar to sco modify for ice shelf
117      END SELECT
118      !
119      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
120         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
121         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
122         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt, Kmm )
123         DEALLOCATE( ztrdu , ztrdv )
124      ENDIF
125      !
126      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg  - Ua: ', mask1=umask,   &
127         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
128      !
129      IF( ln_timing )   CALL timing_stop('dyn_hpg')
130      !
131   END SUBROUTINE dyn_hpg
132
133
134   SUBROUTINE dyn_hpg_init( Kmm )
135      !!----------------------------------------------------------------------
136      !!                 ***  ROUTINE dyn_hpg_init  ***
137      !!
138      !! ** Purpose :   initializations for the hydrostatic pressure gradient
139      !!              computation and consistency control
140      !!
141      !! ** Action  :   Read the namelist namdyn_hpg and check the consistency
142      !!      with the type of vertical coordinate used (zco, zps, sco)
143      !!----------------------------------------------------------------------
144      INTEGER, INTENT( in ) :: Kmm   ! ocean time level index
145      !
146      INTEGER ::   ioptio = 0      ! temporary integer
147      INTEGER ::   ios             ! Local integer output status for namelist read
148      !!
149      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
150      REAL(wp) ::   znad
151      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zts_top, zrhd   ! hypothesys on isf density
152      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  zrhdtop_isf    ! density at bottom of ISF
153      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  ziceload       ! density at bottom of ISF
154      !!
155      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
156         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf
157      !!----------------------------------------------------------------------
158      !
159      READ  ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901)
160901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' )
161      !
162      READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
163902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' )
164      IF(lwm) WRITE ( numond, namdyn_hpg )
165      !
166      IF(lwp) THEN                   ! Control print
167         WRITE(numout,*)
168         WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
169         WRITE(numout,*) '~~~~~~~~~~~~'
170         WRITE(numout,*) '   Namelist namdyn_hpg : choice of hpg scheme'
171         WRITE(numout,*) '      z-coord. - full steps                             ln_hpg_zco    = ', ln_hpg_zco
172         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps
173         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco
174         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf
175         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
176         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj
177      ENDIF
178      !
179      IF( ln_hpg_djc )   &
180         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method',   &
181         &                 '   currently disabled (bugs under investigation).'        ,   &
182         &                 '   Please select either  ln_hpg_sco or  ln_hpg_prj instead' )
183         !
184      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )          &
185         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   &
186         &                 '   the standard jacobian formulation hpg_sco    or '    ,   &
187         &                 '   the pressure jacobian formulation hpg_prj'             )
188         !
189      IF( ln_hpg_isf ) THEN
190         IF( .NOT. ln_isfcav )   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
191       ELSE
192         IF(       ln_isfcav )   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
193      ENDIF
194      !
195      !                               ! Set nhpg from ln_hpg_... flags & consistency check
196      nhpg   = np_ERROR
197      ioptio = 0
198      IF( ln_hpg_zco ) THEN   ;   nhpg = np_zco   ;   ioptio = ioptio +1   ;   ENDIF
199      IF( ln_hpg_zps ) THEN   ;   nhpg = np_zps   ;   ioptio = ioptio +1   ;   ENDIF
200      IF( ln_hpg_sco ) THEN   ;   nhpg = np_sco   ;   ioptio = ioptio +1   ;   ENDIF
201      IF( ln_hpg_djc ) THEN   ;   nhpg = np_djc   ;   ioptio = ioptio +1   ;   ENDIF
202      IF( ln_hpg_prj ) THEN   ;   nhpg = np_prj   ;   ioptio = ioptio +1   ;   ENDIF
203      IF( ln_hpg_isf ) THEN   ;   nhpg = np_isf   ;   ioptio = ioptio +1   ;   ENDIF
204      !
205      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
206      !
207      IF(lwp) THEN
208         WRITE(numout,*)
209         SELECT CASE( nhpg )
210         CASE( np_zco )   ;   WRITE(numout,*) '   ==>>>   z-coord. - full steps '
211         CASE( np_zps )   ;   WRITE(numout,*) '   ==>>>   z-coord. - partial steps (interpolation)'
212         CASE( np_sco )   ;   WRITE(numout,*) '   ==>>>   s-coord. (standard jacobian formulation)'
213         CASE( np_djc )   ;   WRITE(numout,*) '   ==>>>   s-coord. (Density Jacobian: Cubic polynomial)'
214         CASE( np_prj )   ;   WRITE(numout,*) '   ==>>>   s-coord. (Pressure Jacobian: Cubic polynomial)'
215         CASE( np_isf )   ;   WRITE(numout,*) '   ==>>>   s-coord. (standard jacobian formulation) for isf'
216         END SELECT
217         WRITE(numout,*)
218      ENDIF
219      !                         
220   END SUBROUTINE dyn_hpg_init
221
222
223   SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs )
224      !!---------------------------------------------------------------------
225      !!                  ***  ROUTINE hpg_zco  ***
226      !!
227      !! ** Method  :   z-coordinate case, levels are horizontal surfaces.
228      !!      The now hydrostatic pressure gradient at a given level, jk,
229      !!      is computed by taking the vertical integral of the in-situ
230      !!      density gradient along the model level from the suface to that
231      !!      level:    zhpi = grav .....
232      !!                zhpj = grav .....
233      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
234      !!            puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
235      !!            pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
236      !!
237      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
238      !!----------------------------------------------------------------------
239      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
240      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
241      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
242      !
243      INTEGER  ::   ji, jj, jk       ! dummy loop indices
244      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
245      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
246      !!----------------------------------------------------------------------
247      !
248      IF( kt == nit000 ) THEN
249         IF(lwp) WRITE(numout,*)
250         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
251         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
252      ENDIF
253
254      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
255
256      ! Surface value
257      DO_2D_00_00
258         zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
259         ! hydrostatic pressure gradient
260         zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
261         zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
262         ! add to the general momentum trend
263         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
264         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
265      END_2D
266
267      !
268      ! interior value (2=<jk=<jpkm1)
269      DO_3D_00_00( 2, jpkm1 )
270         zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
271         ! hydrostatic pressure gradient
272         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
273            &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
274            &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
275
276         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
277            &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
278            &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
279         ! add to the general momentum trend
280         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
281         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
282      END_3D
283      !
284   END SUBROUTINE hpg_zco
285
286
287   SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs )
288      !!---------------------------------------------------------------------
289      !!                 ***  ROUTINE hpg_zps  ***
290      !!
291      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
292      !!
293      !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
294      !!----------------------------------------------------------------------
295      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
296      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
297      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
298      !!
299      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
300      INTEGER  ::   iku, ikv                         ! temporary integers
301      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
302      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
303      REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv
304      !!----------------------------------------------------------------------
305      !
306      IF( kt == nit000 ) THEN
307         IF(lwp) WRITE(numout,*)
308         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
309         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
310      ENDIF
311
312      ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level
313      CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv )
314
315      ! Local constant initialization
316      zcoef0 = - grav * 0.5_wp
317
318      !  Surface value (also valid in partial step case)
319      DO_2D_00_00
320         zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
321         ! hydrostatic pressure gradient
322         zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
323         zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
324         ! add to the general momentum trend
325         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
326         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
327      END_2D
328
329      ! interior value (2=<jk=<jpkm1)
330      DO_3D_00_00( 2, jpkm1 )
331         zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
332         ! hydrostatic pressure gradient
333         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
334            &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
335            &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
336
337         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
338            &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
339            &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
340         ! add to the general momentum trend
341         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
342         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
343      END_3D
344
345      ! partial steps correction at the last level  (use zgru & zgrv computed in zpshde.F90)
346      DO_2D_00_00
347         iku = mbku(ji,jj)
348         ikv = mbkv(ji,jj)
349         zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj  ,iku,Kmm) )
350         zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji  ,jj+1,ikv,Kmm) )
351         IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
352            puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku)         ! subtract old value
353            zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
354               &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj)
355            puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
356         ENDIF
357         IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
358            pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv)         ! subtract old value
359            zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
360               &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj)
361            pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
362         ENDIF
363      END_2D
364      !
365   END SUBROUTINE hpg_zps
366
367
368   SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs )
369      !!---------------------------------------------------------------------
370      !!                  ***  ROUTINE hpg_sco  ***
371      !!
372      !! ** Method  :   s-coordinate case. Jacobian scheme.
373      !!      The now hydrostatic pressure gradient at a given level, jk,
374      !!      is computed by taking the vertical integral of the in-situ
375      !!      density gradient along the model level from the suface to that
376      !!      level. s-coordinates (ln_sco): a corrective term is added
377      !!      to the horizontal pressure gradient :
378      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
379      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
380      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
381      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
382      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
383      !!
384      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
385      !!----------------------------------------------------------------------
386      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
387      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
388      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
389      !!
390      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
391      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
392      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
393      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj
394      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
395      !!----------------------------------------------------------------------
396      !
397      IF( ln_wd_il ) ALLOCATE(zcpx(jpi,jpj), zcpy(jpi,jpj))
398      !
399      IF( kt == nit000 ) THEN
400         IF(lwp) WRITE(numout,*)
401         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
402         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
403      ENDIF
404      !
405      zcoef0 = - grav * 0.5_wp
406      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
407      ELSE                    ;   znad = 1._wp         ! Variable volume: density
408      ENDIF
409      !
410      IF( ln_wd_il ) THEN
411        DO_2D_00_00
412          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)               ,  ssh(ji+1,jj,Kmm) ) >                &
413               &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
414               &    MAX(  ssh(ji,jj,Kmm) +  ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
415               &                                                       > rn_wdmin1 + rn_wdmin2
416          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)              -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (       &
417               &    MAX(   ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
418               &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
419
420          IF(ll_tmp1) THEN
421            zcpx(ji,jj) = 1.0_wp
422          ELSE IF(ll_tmp2) THEN
423            ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
424            zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
425                        &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) )
426          ELSE
427            zcpx(ji,jj) = 0._wp
428          END IF
429   
430          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
431               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
432               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
433               &                                                      > rn_wdmin1 + rn_wdmin2
434          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        &
435               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
436               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
437
438          IF(ll_tmp1) THEN
439            zcpy(ji,jj) = 1.0_wp
440          ELSE IF(ll_tmp2) THEN
441            ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
442            zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
443                        &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
444          ELSE
445            zcpy(ji,jj) = 0._wp
446          END IF
447        END_2D
448        CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. )
449      END IF
450
451      ! Surface value
452      DO_2D_00_00
453         ! hydrostatic pressure gradient along s-surfaces
454         zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    &
455            &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
456         zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    &
457            &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
458         ! s-coordinate pressure gradient correction
459         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
460            &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
461         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
462            &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
463         !
464         IF( ln_wd_il ) THEN
465            zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
466            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
467            zuap = zuap * zcpx(ji,jj)
468            zvap = zvap * zcpy(ji,jj)
469         ENDIF
470         !
471         ! add to the general momentum trend
472         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap
473         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap
474      END_2D
475
476      ! interior value (2=<jk=<jpkm1)
477      DO_3D_00_00( 2, jpkm1 )
478         ! hydrostatic pressure gradient along s-surfaces
479         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
480            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
481            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
482         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
483            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
484            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
485         ! s-coordinate pressure gradient correction
486         zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
487            &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj)
488         zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
489            &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj)
490         !
491         IF( ln_wd_il ) THEN
492            zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
493            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
494            zuap = zuap * zcpx(ji,jj)
495            zvap = zvap * zcpy(ji,jj)
496         ENDIF
497         !
498         ! add to the general momentum trend
499         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap
500         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap
501      END_3D
502      !
503      IF( ln_wd_il )  DEALLOCATE( zcpx , zcpy )
504      !
505   END SUBROUTINE hpg_sco
506
507
508   SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs )
509      !!---------------------------------------------------------------------
510      !!                  ***  ROUTINE hpg_isf  ***
511      !!
512      !! ** Method  :   s-coordinate case. Jacobian scheme.
513      !!      The now hydrostatic pressure gradient at a given level, jk,
514      !!      is computed by taking the vertical integral of the in-situ
515      !!      density gradient along the model level from the suface to that
516      !!      level. s-coordinates (ln_sco): a corrective term is added
517      !!      to the horizontal pressure gradient :
518      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
519      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
520      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
521      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
522      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
523      !!      iceload is added
524      !!     
525      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
526      !!----------------------------------------------------------------------
527      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
528      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
529      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
530      !!
531      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
532      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
533      REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
534      REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
535      REAL(wp), DIMENSION(jpi,jpj)      ::  zrhdtop_oce
536      !!----------------------------------------------------------------------
537      !
538      zcoef0 = - grav * 0.5_wp   ! Local constant initialization
539      !
540      znad=1._wp                 ! To use density and not density anomaly
541      !
542      !                          ! iniitialised to 0. zhpi zhpi
543      zhpi(:,:,:) = 0._wp   ;   zhpj(:,:,:) = 0._wp
544
545      ! compute rhd at the ice/oce interface (ocean side)
546      ! usefull to reduce residual current in the test case ISOMIP with no melting
547      DO ji = 1, jpi
548        DO jj = 1, jpj
549          ikt = mikt(ji,jj)
550          zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
551          zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
552        END DO
553      END DO
554      CALL eos( zts_top, risfdep, zrhdtop_oce )
555
556!==================================================================================     
557!===== Compute surface value =====================================================
558!==================================================================================
559      DO_2D_00_00
560         ikt    = mikt(ji,jj)
561         iktp1i = mikt(ji+1,jj)
562         iktp1j = mikt(ji,jj+1)
563         ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
564         ! we assume ISF is in isostatic equilibrium
565         zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm)                                    &
566            &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
567            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
568            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
569            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            ) 
570         zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    &
571            &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
572            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
573            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
574            &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            ) 
575         ! s-coordinate pressure gradient correction (=0 if z coordinate)
576         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
577            &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
578         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
579            &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
580         ! add to the general momentum trend
581         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
582         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
583      END_2D
584!==================================================================================     
585!===== Compute interior value =====================================================
586!==================================================================================
587      ! interior value (2=<jk=<jpkm1)
588      DO_3D_00_00( 2, jpkm1 )
589         ! hydrostatic pressure gradient along s-surfaces
590         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
591            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
592            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
593         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
594            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
595            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
596         ! s-coordinate pressure gradient correction
597         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
598            &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj)
599         zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
600            &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj)
601         ! add to the general momentum trend
602         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
603         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
604      END_3D
605      !
606   END SUBROUTINE hpg_isf
607
608
609   SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs )
610      !!---------------------------------------------------------------------
611      !!                  ***  ROUTINE hpg_djc  ***
612      !!
613      !! ** Method  :   Density Jacobian with Cubic polynomial scheme
614      !!
615      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
616      !!----------------------------------------------------------------------
617      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
618      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
619      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
620      !!
621      INTEGER  ::   ji, jj, jk          ! dummy loop indices
622      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
623      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
624      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
625      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
626      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
627      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw
628      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow
629      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k
630      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
631      !!----------------------------------------------------------------------
632      !
633      IF( ln_wd_il ) THEN
634         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
635        DO_2D_00_00
636          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
637               &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
638               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
639               &                                                      > rn_wdmin1 + rn_wdmin2
640          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (        &
641               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                &
642               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
643          IF(ll_tmp1) THEN
644            zcpx(ji,jj) = 1.0_wp
645          ELSE IF(ll_tmp2) THEN
646            ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
647            zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
648                        &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) )
649          ELSE
650            zcpx(ji,jj) = 0._wp
651          END IF
652   
653          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
654               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
655               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
656               &                                                      > rn_wdmin1 + rn_wdmin2
657          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        &
658               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
659               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
660
661          IF(ll_tmp1) THEN
662            zcpy(ji,jj) = 1.0_wp
663          ELSE IF(ll_tmp2) THEN
664            ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
665            zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
666                        &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
667          ELSE
668            zcpy(ji,jj) = 0._wp
669          END IF
670        END_2D
671        CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. )
672      END IF
673
674      IF( kt == nit000 ) THEN
675         IF(lwp) WRITE(numout,*)
676         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
677         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
678      ENDIF
679
680      ! Local constant initialization
681      zcoef0 = - grav * 0.5_wp
682      z1_10  = 1._wp / 10._wp
683      z1_12  = 1._wp / 12._wp
684
685      !----------------------------------------------------------------------------------------
686      !  compute and store in provisional arrays elementary vertical and horizontal differences
687      !----------------------------------------------------------------------------------------
688
689!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
690
691      DO_3D_00_00( 2, jpkm1 )
692         drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
693         dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1)
694         drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
695         dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  )
696         drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
697         dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  )
698      END_3D
699
700      !-------------------------------------------------------------------------
701      ! compute harmonic averages using eq. 5.18
702      !-------------------------------------------------------------------------
703      zep = 1.e-15
704
705!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
706!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
707
708      DO_3D_00_00( 2, jpkm1 )
709         cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
710
711         cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  )
712         cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  )
713
714         cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  )
715         cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  )
716
717         IF( cffw > zep) THEN
718            drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   &
719               &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
720         ELSE
721            drhow(ji,jj,jk) = 0._wp
722         ENDIF
723
724         dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   &
725            &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
726
727         IF( cffu > zep ) THEN
728            drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   &
729               &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
730         ELSE
731            drhou(ji,jj,jk ) = 0._wp
732         ENDIF
733
734         IF( cffx > zep ) THEN
735            dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   &
736               &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
737         ELSE
738            dzu(ji,jj,jk) = 0._wp
739         ENDIF
740
741         IF( cffv > zep ) THEN
742            drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   &
743               &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
744         ELSE
745            drhov(ji,jj,jk) = 0._wp
746         ENDIF
747
748         IF( cffy > zep ) THEN
749            dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   &
750               &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
751         ELSE
752            dzv(ji,jj,jk) = 0._wp
753         ENDIF
754
755      END_3D
756
757      !----------------------------------------------------------------------------------
758      ! apply boundary conditions at top and bottom using 5.36-5.37
759      !----------------------------------------------------------------------------------
760      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  )
761      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  )
762      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  )
763
764      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
765      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
766      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
767
768
769      !--------------------------------------------------------------
770      ! Upper half of top-most grid box, compute and store
771      !-------------------------------------------------------------
772
773!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified
774!          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
775
776      DO_2D_00_00
777         rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               &
778            &                   * (  rhd(ji,jj,1)                                     &
779            &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  &
780            &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  &
781            &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  )
782      END_2D
783
784!!bug gm    : here also, simplification is possible
785!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
786
787      DO_3D_00_00( 2, jpkm1 )
788
789         rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
790            &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   &
791            &            - grav * z1_10 * (                                                                   &
792            &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     &
793            &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
794            &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     &
795            &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   &
796            &                             )
797
798         rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   &
799            &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   &
800            &            - grav* z1_10 * (                                                                    &
801            &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     &
802            &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
803            &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     &
804            &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   &
805            &                            )
806
807         rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 &
808            &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   &
809            &            - grav* z1_10 * (                                                                    &
810            &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     &
811            &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
812            &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     &
813            &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   &
814            &                            )
815
816      END_3D
817      CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1., rho_i, 'U', 1., rho_j, 'V', 1. )
818
819      ! ---------------
820      !  Surface value
821      ! ---------------
822      DO_2D_00_00
823         zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
824         zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
825         IF( ln_wd_il ) THEN
826           zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
827           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
828         ENDIF
829         ! add to the general momentum trend
830         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
831         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
832      END_2D
833
834      ! ----------------
835      !  interior value   (2=<jk=<jpkm1)
836      ! ----------------
837      DO_3D_00_00( 2, jpkm1 )
838         ! hydrostatic pressure gradient along s-surfaces
839         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
840            &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
841            &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj)
842         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
843            &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
844            &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
845         IF( ln_wd_il ) THEN
846           zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
847           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
848         ENDIF
849         ! add to the general momentum trend
850         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
851         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
852      END_3D
853      !
854      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
855      !
856   END SUBROUTINE hpg_djc
857
858
859   SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs )
860      !!---------------------------------------------------------------------
861      !!                  ***  ROUTINE hpg_prj  ***
862      !!
863      !! ** Method  :   s-coordinate case.
864      !!      A Pressure-Jacobian horizontal pressure gradient method
865      !!      based on the constrained cubic-spline interpolation for
866      !!      all vertical coordinate systems
867      !!
868      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
869      !!----------------------------------------------------------------------
870      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
871      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
872      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
873      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
874      !!
875      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
876      REAL(wp) ::   zcoef0, znad                    ! local scalars
877      !
878      !! The local variables for the correction term
879      INTEGER  :: jk1, jis, jid, jjs, jjd
880      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables
881      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
882      REAL(wp) :: zrhdt1
883      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
884      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
885      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
886      REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n
887      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
888      !!----------------------------------------------------------------------
889      !
890      IF( kt == nit000 ) THEN
891         IF(lwp) WRITE(numout,*)
892         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
893         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
894      ENDIF
895
896      ! Local constant initialization
897      zcoef0 = - grav
898      znad = 1._wp
899      IF( ln_linssh )   znad = 0._wp
900
901      IF( ln_wd_il ) THEN
902         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
903         DO_2D_00_00
904          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
905               &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
906               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
907               &                                                      > rn_wdmin1 + rn_wdmin2
908          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (         &
909               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                &
910               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
911
912          IF(ll_tmp1) THEN
913            zcpx(ji,jj) = 1.0_wp
914          ELSE IF(ll_tmp2) THEN
915            ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
916            zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
917                        &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
918           
919             zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
920          ELSE
921            zcpx(ji,jj) = 0._wp
922          END IF
923   
924          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
925               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
926               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
927               &                                                      > rn_wdmin1 + rn_wdmin2
928          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (      &
929               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
930               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
931
932          IF(ll_tmp1) THEN
933            zcpy(ji,jj) = 1.0_wp
934          ELSE IF(ll_tmp2) THEN
935            ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
936            zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
937                        &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
938             zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
939
940            ELSE
941               zcpy(ji,jj) = 0._wp
942            ENDIF
943         END_2D
944         CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1., zcpy, 'V', 1. )
945      ENDIF
946
947      ! Clean 3-D work arrays
948      zhpi(:,:,:) = 0._wp
949      zrhh(:,:,:) = rhd(:,:,:)
950
951      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
952      DO_2D_11_11
953       jk = mbkt(ji,jj)+1
954       IF(     jk <=  0   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
955       ELSEIF( jk ==  1   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
956       ELSEIF( jk < jpkm1 ) THEN
957          DO jkk = jk+1, jpk
958             zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   &
959                &                      gde3w(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2))
960          END DO
961       ENDIF
962      END_2D
963
964      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
965      DO_2D_11_11
966         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad
967      END_2D
968
969      DO_3D_11_11( 2, jpk )
970         zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm)
971      END_3D
972
973      fsp(:,:,:) = zrhh (:,:,:)
974      xsp(:,:,:) = zdept(:,:,:)
975
976      ! Construct the vertical density profile with the
977      ! constrained cubic spline interpolation
978      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3
979      CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
980
981      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
982      DO_2D_01_01
983       zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
984          &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
985
986       ! assuming linear profile across the top half surface layer
987       zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1
988      END_2D
989
990      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
991      DO_3D_01_01( 2, jpkm1 )
992      zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  &
993         &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   &
994         &                           asp  (ji,jj,jk-1), bsp  (ji,jj,jk-1), &
995         &                           csp  (ji,jj,jk-1), dsp  (ji,jj,jk-1)  )
996      END_3D
997
998      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
999
1000      ! Prepare zsshu_n and zsshv_n
1001      DO_2D_00_00
1002!!gm BUG ?    if it is ssh at u- & v-point then it should be:
1003!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
1004!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
1005!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * &
1006!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
1007!!gm not this:
1008       zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
1009                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
1010       zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
1011                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
1012      END_2D
1013
1014      CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1., zsshv_n, 'V', 1. )
1015
1016      DO_2D_00_00
1017       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 
1018       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)
1019      END_2D
1020
1021      DO_3D_00_00( 2, jpkm1 )
1022      zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
1023      zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
1024      END_3D
1025
1026      DO_3D_00_00( 1, jpkm1 )
1027      zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
1028      zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
1029      END_3D
1030
1031      DO_3D_00_00( 1, jpkm1 )
1032      zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
1033      zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
1034      zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
1035      zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
1036      END_3D
1037
1038
1039      DO_3D_00_00( 1, jpkm1 )
1040      zpwes = 0._wp; zpwed = 0._wp
1041      zpnss = 0._wp; zpnsd = 0._wp
1042      zuijk = zu(ji,jj,jk)
1043      zvijk = zv(ji,jj,jk)
1044
1045      !!!!!     for u equation
1046      IF( jk <= mbku(ji,jj) ) THEN
1047         IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
1048           jis = ji + 1; jid = ji
1049         ELSE
1050           jis = ji;     jid = ji +1
1051         ENDIF
1052
1053         ! integrate the pressure on the shallow side
1054         jk1 = jk
1055         DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
1056           IF( jk1 == mbku(ji,jj) ) THEN
1057             zuijk = -zdept(jis,jj,jk1)
1058             EXIT
1059           ENDIF
1060           zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
1061           zpwes = zpwes +                                    &
1062                integ_spline(zdept(jis,jj,jk1), zdeps,            &
1063                       asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
1064                       csp(jis,jj,jk1),    dsp(jis,jj,jk1))
1065           jk1 = jk1 + 1
1066         END DO
1067
1068         ! integrate the pressure on the deep side
1069         jk1 = jk
1070         DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
1071           IF( jk1 == 1 ) THEN
1072             zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad)
1073             zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
1074                                               bsp(jid,jj,1),   csp(jid,jj,1), &
1075                                               dsp(jid,jj,1)) * zdeps
1076             zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
1077             EXIT
1078           ENDIF
1079           zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
1080           zpwed = zpwed +                                        &
1081                  integ_spline(zdeps,              zdept(jid,jj,jk1), &
1082                         asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
1083                         csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
1084           jk1 = jk1 - 1
1085         END DO
1086
1087         ! update the momentum trends in u direction
1088
1089         zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
1090         IF( .NOT.ln_linssh ) THEN
1091           zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
1092              &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) )
1093          ELSE
1094           zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
1095         ENDIF
1096         IF( ln_wd_il ) THEN
1097            zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)
1098            zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
1099         ENDIF
1100         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
1101      ENDIF
1102
1103      !!!!!     for v equation
1104      IF( jk <= mbkv(ji,jj) ) THEN
1105         IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
1106           jjs = jj + 1; jjd = jj
1107         ELSE
1108           jjs = jj    ; jjd = jj + 1
1109         ENDIF
1110
1111         ! integrate the pressure on the shallow side
1112         jk1 = jk
1113         DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
1114           IF( jk1 == mbkv(ji,jj) ) THEN
1115             zvijk = -zdept(ji,jjs,jk1)
1116             EXIT
1117           ENDIF
1118           zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
1119           zpnss = zpnss +                                      &
1120                  integ_spline(zdept(ji,jjs,jk1), zdeps,            &
1121                         asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
1122                         csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
1123           jk1 = jk1 + 1
1124         END DO
1125
1126         ! integrate the pressure on the deep side
1127         jk1 = jk
1128         DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
1129           IF( jk1 == 1 ) THEN
1130             zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)
1131             zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
1132                                               bsp(ji,jjd,1),   csp(ji,jjd,1), &
1133                                               dsp(ji,jjd,1) ) * zdeps
1134             zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
1135             EXIT
1136           ENDIF
1137           zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
1138           zpnsd = zpnsd +                                        &
1139                  integ_spline(zdeps,              zdept(ji,jjd,jk1), &
1140                         asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1141                         csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1142           jk1 = jk1 - 1
1143         END DO
1144
1145
1146         ! update the momentum trends in v direction
1147
1148         zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
1149         IF( .NOT.ln_linssh ) THEN
1150            zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
1151                    ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )
1152         ELSE
1153            zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
1154         ENDIF
1155         IF( ln_wd_il ) THEN
1156            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 
1157            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 
1158         ENDIF
1159
1160         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
1161      ENDIF
1162         !
1163      END_3D
1164      !
1165      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
1166      !
1167   END SUBROUTINE hpg_prj
1168
1169
1170   SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
1171      !!----------------------------------------------------------------------
1172      !!                 ***  ROUTINE cspline  ***
1173      !!
1174      !! ** Purpose :   constrained cubic spline interpolation
1175      !!
1176      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3
1177      !!
1178      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
1179      !!----------------------------------------------------------------------
1180      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   fsp, xsp           ! value and coordinate
1181      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   asp, bsp, csp, dsp ! coefficients of the interpoated function
1182      INTEGER                   , INTENT(in   ) ::   polynomial_type    ! 1: cubic spline   ;   2: Linear
1183      !
1184      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
1185      INTEGER  ::   jpi, jpj, jpkm1
1186      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
1187      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha
1188      REAL(wp) ::   zdf(size(fsp,3))
1189      !!----------------------------------------------------------------------
1190      !
1191!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!! 
1192      jpi   = size(fsp,1)
1193      jpj   = size(fsp,2)
1194      jpkm1 = MAX( 1, size(fsp,3) - 1 )
1195      !
1196      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
1197         DO ji = 1, jpi
1198            DO jj = 1, jpj
1199           !!Fritsch&Butland's method, 1984 (preferred, but more computation)
1200           !    DO jk = 2, jpkm1-1
1201           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)
1202           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1203           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
1204           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
1205           !
1206           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
1207           !
1208           !       IF(zdf1 * zdf2 <= 0._wp) THEN
1209           !           zdf(jk) = 0._wp
1210           !       ELSE
1211           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
1212           !       ENDIF
1213           !    END DO
1214
1215           !!Simply geometric average
1216               DO jk = 2, jpkm1-1
1217                  zdf1 = (fsp(ji,jj,jk  ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk  ) - xsp(ji,jj,jk-1))
1218                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk  )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk  ))
1219
1220                  IF(zdf1 * zdf2 <= 0._wp) THEN
1221                     zdf(jk) = 0._wp
1222                  ELSE
1223                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
1224                  ENDIF
1225               END DO
1226
1227               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
1228                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) )           -  0.5_wp * zdf(2)
1229               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
1230                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1)
1231
1232               DO jk = 1, jpkm1 - 1
1233                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1234                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1235                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1236                 zddf1  = -2._wp * ztmp1 + ztmp2
1237                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1238                 zddf2  =  2._wp * ztmp1 - ztmp2
1239
1240                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1241                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1242                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1243                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1244                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
1245                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
1246                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
1247                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
1248                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
1249               END DO
1250            END DO
1251         END DO
1252
1253      ELSEIF ( polynomial_type == 2 ) THEN     ! Linear
1254         DO ji = 1, jpi
1255            DO jj = 1, jpj
1256               DO jk = 1, jpkm1-1
1257                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1258                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1259
1260                  dsp(ji,jj,jk) = 0._wp
1261                  csp(ji,jj,jk) = 0._wp
1262                  bsp(ji,jj,jk) = ztmp1 / zdxtmp
1263                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1264               END DO
1265            END DO
1266         END DO
1267         !
1268      ELSE
1269         CALL ctl_stop( 'invalid polynomial type in cspline' )
1270      ENDIF
1271      !
1272   END SUBROUTINE cspline
1273
1274
1275   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)
1276      !!----------------------------------------------------------------------
1277      !!                 ***  ROUTINE interp1  ***
1278      !!
1279      !! ** Purpose :   1-d linear interpolation
1280      !!
1281      !! ** Method  :   interpolation is straight forward
1282      !!                extrapolation is also permitted (no value limit)
1283      !!----------------------------------------------------------------------
1284      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr
1285      REAL(wp)             ::  f ! result of the interpolation (extrapolation)
1286      REAL(wp)             ::  zdeltx
1287      !!----------------------------------------------------------------------
1288      !
1289      zdeltx = xr - xl
1290      IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN
1291         f = 0.5_wp * (fl + fr)
1292      ELSE
1293         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1294      ENDIF
1295      !
1296   END FUNCTION interp1
1297
1298
1299   FUNCTION interp2( x, a, b, c, d )  RESULT(f)
1300      !!----------------------------------------------------------------------
1301      !!                 ***  ROUTINE interp1  ***
1302      !!
1303      !! ** Purpose :   1-d constrained cubic spline interpolation
1304      !!
1305      !! ** Method  :  cubic spline interpolation
1306      !!
1307      !!----------------------------------------------------------------------
1308      REAL(wp), INTENT(in) ::  x, a, b, c, d
1309      REAL(wp)             ::  f ! value from the interpolation
1310      !!----------------------------------------------------------------------
1311      !
1312      f = a + x* ( b + x * ( c + d * x ) )
1313      !
1314   END FUNCTION interp2
1315
1316
1317   FUNCTION interp3( x, a, b, c, d )  RESULT(f)
1318      !!----------------------------------------------------------------------
1319      !!                 ***  ROUTINE interp1  ***
1320      !!
1321      !! ** Purpose :   Calculate the first order of derivative of
1322      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3
1323      !!
1324      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2
1325      !!
1326      !!----------------------------------------------------------------------
1327      REAL(wp), INTENT(in) ::  x, a, b, c, d
1328      REAL(wp)             ::  f ! value from the interpolation
1329      !!----------------------------------------------------------------------
1330      !
1331      f = b + x * ( 2._wp * c + 3._wp * d * x)
1332      !
1333   END FUNCTION interp3
1334
1335
1336   FUNCTION integ_spline( xl, xr, a, b, c, d )  RESULT(f)
1337      !!----------------------------------------------------------------------
1338      !!                 ***  ROUTINE interp1  ***
1339      !!
1340      !! ** Purpose :   1-d constrained cubic spline integration
1341      !!
1342      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr
1343      !!
1344      !!----------------------------------------------------------------------
1345      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d
1346      REAL(wp)             ::  za1, za2, za3
1347      REAL(wp)             ::  f                   ! integration result
1348      !!----------------------------------------------------------------------
1349      !
1350      za1 = 0.5_wp * b
1351      za2 = c / 3.0_wp
1352      za3 = 0.25_wp * d
1353      !
1354      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1355         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1356      !
1357   END FUNCTION integ_spline
1358
1359   !!======================================================================
1360END MODULE dynhpg
1361
Note: See TracBrowser for help on using the repository browser.