New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dynhpg.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DYN/dynhpg.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

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