source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynhpg.F90 @ 12546

Last change on this file since 12546 was 12546, checked in by orioltp, 6 months ago

Adding precision specification in hardcoded reals and other modifications to allow compilation without forcing reals without precision specification to a certain value through compiler flags

  • 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.0_wp, zcpy, 'V', 1.0_wp )
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.0_wp, zcpy, 'V', 1.0_wp )
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.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp )
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.0_wp, zcpy, 'V', 1.0_wp )
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.0_wp, zsshv_n, 'V', 1.0_wp )
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.