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/2020/dev_r13787_doc_latex_recovery/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/DYN/dynhpg.F90 @ 14066

Last change on this file since 14066 was 14066, checked in by nicolasmartin, 3 years ago
#2414 Sync merge with trunk
Property svn:keywords set to `Id`
File size: 75.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: